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Abstract 



The majority of galactic gamma rays are produced by interaction of cosmic rays 
with matter. This results in a diffuse radiation concentrated in the galactic plane 
where the flux of cosmic rays and the density of material (mostly atomic, molecular 
and ionized hydrogen) is high. The interactions producing gamma rays include, 
among others, the decay of 7r°'s produced in spallation reactions. Gamma emission 
from the plane has indeed been detected in the energy range up to 30 GeV by 
space-based detectors. Above 1 GeV, the observed intensity is notably higher 
than expected in simple models, possibly implying an enhancement at the TeV 
region as well. Observations at TeV energies, for which the flux is too low for 
satellite detection, can be done with ground based telescopes. Milagro is a large 
aperture water Cherenkov detector for extensive air showers, collecting data from 
a solid angle of more than two steradians in the overhead sky at energies near 1 
TeV. A 2000-2001 data set from Milagro has been used to search for the emission 
of diffuse gamma rays from the galactic disk. An excess has been observed from 
the region of the Milagro inner Galaxy defined by I G (20°, 100°) and |6| < 5° 
with the significance 2.3 ■ 10 -4 . The emission from the region of the Milagro outer 
Galaxy defined by / € (140°, 220°) and |6| < 5° is not inconsistent with being 
that of background only. Under the assumption that EGRET measurements in 
10-30 GeV range can be extended to TeV region with a simple power law energy 
spectrum, the integral gamma ray flux with energies above 1 TeV for the region of 
inner Galaxy is measured to be F(> ITeV) = (9.5±2.0)-10~ 10 cmT x sr' 1 s" 1 with 
spectral index a 1 = 2.59 ± 0.07. The 99.9% upper limit for the diffuse emission in 
the region of outer Galaxy is set at F(> ITeV) < 4.5 • 10~ 10 cm -1 sr _1 s" 1 using 
a differential spectral index of a 7 = 2.49. The upper limit for the outer Galaxy 
is consistent with the extrapolation of EGRET measurements between 1 and 30 
GeV. Extrapolation of the EGRET measurements between 1 and 30 GeV for the 
region of inner Galaxy using constant power law spectral index is incompatible 
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with the Milagro data. This indicates softening of the spectrum at energy between 
10 GeV and 1 TeV. These observations may be used to constrain some models of 
Galactic gamma ray emission. 
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Chapter 1 
Introduction 



Most cosmic rays are accelerated by unknown objects in our Galaxy and are 
trapped (for about 100 million years) by Galactic magnetic fields. The interac- 
tion of high energy cosmic rays with the interstellar material produces 7-rays by 
a combination of electron bremsstrahlung, inverse Compton and nucleon-nucleon 
processes. The nucleon-nucleon interactions give rise to 7r°'s which decay to gamma 
rays and are expected to dominate the flux at energies above several GeV. In this 
manner, the regions of enhanced density (clouds of mostly atomic and molecular 
hydrogen) act as passive targets, converting some fraction of impinging cosmic 
rays into gamma rays. This should appear as a diffuse glow concentrated in the 
narrow band along the Galactic equator. Indeed, such an emission was detected 
by the space-borne detectors SAS 2 [6], COS B [5] and EGRET [7] at energies 
up to 30 GeV. Figure 1.1 presents the EGRET all sky survey plotted in Galactic 
coordinates. The Galactic plane is clearly visible. 



However, observations with present satellite based instruments at higher ener- 
gies are not possible due to the rapidly decreasing flux of 7-rays, requiring bigger 
effective area of the detectors. Therefore, the use of ground-based arrays is need- 
ed to observe the diffuse Galactic radiation. Inasmuch as the Galactic cosmic ray 
spectrum extends beyond 10 15 eV, the diffuse Galactic emission should extend well 
beyond the energy threshold of Milagro (~ 400 GeV). A number of authors have 
estimated the expected diffuse very high energy gamma-ray flux from the Galactic 
plane (see for example [4]): they generally predict a flux within ±5° of the Galactic 
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EGRET All-Sky Gamma-Ray Survey Above 100 MeV 




Figure 1.1: All sky gamma ray survey presented in Galactic coordinates. 
Galactic center is in the middle of the map [10]. 



equator in latitude that is ~ 10~ 4 — 10~ 5 of the cosmic ray flux for the regions 
of the outer galaxy. 1 The shape of the gamma-ray spectrum is predicted by the 
same authors to have power law form ^ ~ E~ a , with spectral index a = 2.7. 
However, at TeV energies the contribution from source cosmic rays, considered by 
[2], may increase the expected diffuse 7-ray flux by almost an order of magnitude 
compared to 7r°-decay model predictions. It is also possible that the spectrum of 
cosmic rays in the interstellar medium is substantially harder compared with the 
local one measured directly in the solar neighborhood [1] which will lead to higher 
diffuse 7-ray flux as well. 

At present, gamma rays from the galactic plane have not been detected above 
EGRET energies (only upper limits were set). The only measurements that ap- 
proach the required sensitivity are above 180 TeV, performed by the CASA-MIA 
experiment [3]. The best measurement in the 1 TeV region, which is two orders 
of magnitude less sensitive, is due to Whipple [8]. The present state of theoretical 
predictions and experimental measurements is summarized in figure 1.2 [9]. Mila- 
gro, a detector designed to cover the energy gap in the few TeV region between 
other existing instruments, should be able to detect the diffuse very high energy 
Galactic emission and possibly its spatial distribution and provide an enhanced 
understanding of Galactic cosmic rays. The sky coverage of Milagro is illustrated 
in figure 1.3. Because Milagro is located in the northern hemisphere at a latitude of 
36°, the Galactic center is not in its field of view. However, a considerable portion 
of the outer disk is visible to Milagro. For a year's exposure, Milagro is sensitive 
to a gamma ray flux of about that theoretically predicted. 



1 The outer Galaxy is denned as the region with galactic longitude I, 40° < I < 320°. 



3 



10" 



10" 4 



Galactic plane 



10 



-7 



-8 



10 



10"' 




- Outer Galaxy 



10 



10 



3 10 



-11 



10" 



12 



c 10 



-13 



10 



10 



10 



-14 



15 



16 



ii] — I i rmn| — i i i iiuij — i i nni i[ — i i hu b 

A Whipple 
O Tibet 

V HEGRA-AIROBICC 
□ EAS-TOP 
* Utah-Mich 
O CASA-MIA 



t 



i i mill — i i i mill i i i mill i i i mill i i i mill 



T 



■t 



I I I llllll ' i i mill 



10 1 10 2 10 3 10 4 
Energy (GeV) 



10 s 



10° 



10 7 



Figure 1.2: Theoretical predictions and observations of the diffuse 7-ray flux from 
the Galactic plane as a function of energy. EGRET measurements for the inner and 
outer regions are shown by filled points with error bars. The solid line represents 
the fit to the EGRET data for the inner Galaxy indicating that the observed 
spectrum is significantly harder than expected (dashed line). The dotted line 
shows predictions for the outer Galaxy region. Identification of other instruments 
and their upper limits are shown as open symbols [9]. 
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Figure 1.3: The density of events from Milagro (in arbitrary units) plotted in a 
galactic coordinate system, sine equal area projection. Grid lines are plotted every 
30° in longitude and latitude. The Galactic center, not visible by Milagro, is in 
the center of the map. Galactic longitude increases to the left. 
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Chapter 2 



Diffuse Galactic Gamma Ray 
Emission. 

High energy gamma rays produced by interactions of cosmic rays with the inter- 
stellar material should provide tracers of cosmic rays because their trajectories are 
undeflected by interstellar magnetic fields and because they traverse the Galaxy 
without significant attenuation. (Neutrinos would provide even better tracers be- 
cause of their weak interactions.) Therefore, diffuse gamma ray emission of the 
Galactic disk carries unique information about the production sites, the fluxes and 
the spatial distribution of Galactic cosmic rays. Indeed, the separation of different 
emission mechanisms in a broad energy region from a few KeV to few hundred 
TeV in different parts of the Galaxy would provide an important insight into the 
problem of the origin of cosmic rays and of their propagation in the interstellar 
medium. This is illustrated on figure 2.1 [1]. Diffuse gamma radiation in the cen- 
ter of the Galaxy or in its halo has been proposed as a probe of annihilating dark 
matter particles as well [11]. 

The observations of the diffuse gamma ray radiation conducted in 1990's by the 
EGRET detector aboard Compton Gamma Ray Observatory [7] resulted in good 
quality data over three decades in energy of gamma rays. The results support 
a Galactic origin of cosmic rays and strong correlation between the high energy 
gamma ray flux and the column density of the galactic hydrogen. The latter 
demonstrated the existence of a truly diffuse radiation based on the earlier data 
from SAS-2 and COS B [12]. 

The extraction of the diffuse component of the gamma ray emission from the 
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EGRET 




Log(E/sV) 

(a) The heavy solid line is the flux with- 
out contribution from positron annihila- 
tion, heavy dashed line takes it into ac- 
count. 




U>g(6/4V) 



(b) The heavy solid line is the flux from 
the main component (figure 2.1(a)), 
heavy dashed line takes also into ac- 
count a possible population of electrons 
accelerated to energies 250 TeV. Mea- 
surements from X-rays to very high en- 
ergy gamma rays and upper limits are 
also shown on the plot. 



Figure 2.1: The fluxes of diffuse radiation produced by both electronic and nucle- 
onic components of cosmic rays in the inner Galaxy (315° < / < 45°) calculated for 
a hard power law spectrum of electrons (with index 2.15) and protons (with index 
2.1 and gradual turnover to 2.75 above few GeV). Contributions from different 
emission mechanisms are shown and labeled. (Adopted from [1]) 
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Figure 2.2: Longitude dependence of diffuse gamma ray emission averaged over 
4° wide latitude intervals in the energy range between 1 and 30 GeV measured 
by EGRET. The emission including point sources is shown as dotted line. The 
solid line is the best fit model calculated for the same energy range, convolved 
with the EGRET resolution function and averaged over the same latitude interval. 
The isotropic diffuse component added to the model for this energy range is 0.12 • 
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The model without the 60 % increase is shown as a 



dashed line. (Adopted from [7].) 
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Figure 2.3: Latitude dependence of diffuse gamma ray emission averaged over 20° 
wide longitude intervals in the energy range between 1 and 30 GeV measured by 
EGRET. See caption of figure 2.2 for details. 
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100 1000 10000 

Energy [MeV] 

Figure 2.4: Average diffuse gamma ray spectrum of inner Galaxy region, 300° < I < 
60°, \b\ < 10° measured by EGRET. The contribution from point sources had been 
removed. The model calculations plus the isotropic diffuse emission (dash-dotted 
line ID) is shown as solid line. The individual components of the calculation are 
also shown as dashed lines: nucleon-nucleaon (NN), electron bremsstrahlung (EB), 
inverse Compton scattering of electrons(IC) contributions. At energies above 1 
GeV the data points are significantly higher than the model calculations. (Adopted 
from [7].) 
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EGRET data, that is the radiation produced by cosmic ray electrons, protons 
and nuclei interacting with the ambient gas and photon fields, is obscured by the 
emission from resolved and an uncertain number of unresolved point sources as 
well as by the diffuse isotropic emission of presumably extragalactic origin. An 
example of longitudinal intensity distribution of the observed emission, including 
the isotropic emission and excluding the point source contribution is shown on 
figure 2.2. Figure 2.3 shows the latitude distribution of the intensity for the same 
energy range. The dotted lines represent the observed emission including the point 
sources, the solid lines are the calculated intensities using the model described in 
[7, 14]. In this model, the EGRET data together with radio data was used to 
develop a three-dimensional picture of both gas and cosmic ray densities in the 
Galaxy. The diffuse gamma ray emission for energy E 1 from the galactic longitude 
/ and latitude b, \b\ < 10° is given by [7, 14]: 

j(E y ,l,b) = J [c e (p,l,b)q em (E 1 ) + c n (p,l,b)q nm (E J )} x 
x [n H i(p, I, b) + n H n(p, I, b) + n H2 (p, I, b)] dp+ 

+ ^E/ c e (p,l,b)qU(E„p)u$\p,l,b)dp 

(photons cm~ 2 sr~ x GeV^ 1 ) 

where p is the distance of the line of sight in the direction of I and b mea- 
sured from the Sun. The first integral represents the gamma ray production 
due to cosmic ray interactions with matter where q em {E^) and q nm (E^) are the 
electron bremsstrahlung and nucleon-nucleon production functions per target hy- 
drogen atom based on the cosmic ray spectrum measured in the vicinity of the 
Sun. The functions c e (p,l,b) and c n (p,l,b) are the ratios of the cosmic ray elec- 
tron and proton densities at the given location to their densities in the vicin- 
ity of the Sun. The electron and proton spectra are assumed to have the same 
shape as measured locally near the Sun and therefore the c's are independent 
of energy. The ratios are also assumed to be be equal and independent of b: 
c e (p,l,b) = c n (p,l,b) = c(p, /), c(p = 0, 1) = 1. The quantities nni(p,l,b) and 
nH 2 (p,l,b) are the atomic and molecular hydrogen densities expressed as atoms 
per unit volume derived from the 21 cm hyperfine transition emission line (HI) 
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and 2.6 mm kinematic J = 1 — > transition line of CO surveys. The CO inten- 
sities are scaled by 2X where X = N(H 2 )/Wco is the proportionality constant 
between the column density of molecular hydrogen and the integrated intensity 
of the CO line. The distribution of ionized hydrogen uhii is taken from a model 
[7] and is shown to have small contribution to the diffuse gamma ray emission 
compared to that of HI and H 2 . The second integral describes the contribution 
from inverse Compton interactions between cosmic ray electrons and interstellar 
photons where q^(E y , p) is the inverse Compton production function based on the 
local electron spectrum, and the summation is over discrete wavelength bands (i) of 
cosmic blackbody radiation, infrared, optical and ultraviolet that arise from within 
our Galaxy with corresponding photon energy density distributions u^\p, I, b). 

This model accurately matches the observed emission as seen by EGRET for 
all Galactic longitudes over the energy range from 30 MeV to about 30 GeV. The 
model underestimates the observations at energies above 1 GeV in the region of 
inner Galaxy 300° < I < 60°, |6| < 10°. (The model calculations without the 
60 % increase in this region is shown as a dashed line on figures 2.2, 2.3. This 
is also illustrated on figure 2.4.) Some authors [13] have suggested that because 
electron propagation is limited by radiative losses the local measurements of the 
electron spectrum may not be representative for the entire Galaxy. In this 
harder cosmic ray electron spectrum could be used to explain the observed excess. 
It is also possible [1] that the proton spectrum in the inner Galaxy is harder then 
observed in the Solar neighborhood, at least in the region below few GeV. Another 
possibility is contribution from unresolved sources. Assuming these are supernova 
remnants, it was shown [2] that their spatially averaged contribution to the diffuse 
gamma ray flux at 1 TeV should exceed the above model predictions [7] extended 
to 1 TeV by almost an order of magnitude. It is therefore of interest to search for 
the diffuse emission from the Galactic plane at energies around 1 TeV. 
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Chapter 3 



Extensive Air Showers. 



3.1 Longitudinal Development of Extensive Air 
Showers. 

The desire to detect low fluxes of gamma rays at energy around 1 TeV using ground 
based telescopes impels consideration of propagation of the photons to the ground 
level, where they can be registered. 

A high energy primary gamma-ray entering the atmosphere interacts with it, 
initiating the production of secondary particles which in turn create tertiary par- 
ticles and so on. Such electromagnetic cascades are propagated predominantly by 
photons and electrons. The basic high energy processes that make up the cascade 
are pair production and bremsstrahlung occurring in the field of nuclei of air which 
produce successive generations of electrons, positrons and photons. Charged par- 
ticles are removed from the shower development by ionization losses, photons - 
by Compton scattering. The number of particles increases until their energies de- 
crease to the critical energy E c 80MeV, when ionization and scattering become 
the main energy loss mechanisms. This stage of the development is called shower 
maximum. 

Inasmuch as for ultrarelativistic particles the radiation length X for the brems- 
strahlung process (the amount of the material traversed over which particle's en- 
ergy is decreased by a factor of e, E — E e~ x / X °) is approximately equal to the 
gamma ray interaction length for electron-positron pair production [15], it provides 
a convenient scale, and in the case of air it is about X 37g/cm 2 . 
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The cascade development provides an interesting example of a stochastic pro- 
cess which is prohibitively difficult for analytic calculations. Although several 
approximations have been considered [16], numerical methods are usually needed 
to obtain results of practical use. 

Recent simulations by DiSciascio et al [17] show that the average number of 
photons and electrons with the energy greater than E t h in the shower initiated by 
a photon with energy E is well described by a modified Greisen formula: 



31 

N(E ,E th ,t) = A(E th )——e tl( - 1 ~ 1 ' 5lns ^ (3.1) 

Vv 

where t\ is the modified depth t = x/X from the top of the atmosphere 
measured along the trajectory of the primary particle and expressed in the units 
of radiation length: 



h = t + a{E th ) 



The parameter s\ represents the age of the shower and increases as it develops 
starting at s± — 0, s± — 1 at the maximum, si > 1 in the declining stage of the 
cascade: 



Sl 



3ti 



h + 2y 



y = in- 



c 



The parameterization is valid in the depth range 4 < t < 24 for primary photon 
energies 0.1 < E < 10 3 TeV. The coefficients A(E th ) and a(E th ) are given in the 
table 3.1. The dependence is illustrated on figure 3.1. The Milagro detector is 
situated at the vertical depth of about 20 radiation lengths, thus within the range 
of the simulations for photons with zenith angles between and 32 degrees. 
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Eth, 
MeV 


electrons 


photons 


A 


a 


A 


a 


1 


0.92 


0.00 


4.80 


-0.88 


5 


0.75 


0.19 


2.98 


-0.69 


10 


0.63 


0.35 


2.13 


-0.57 


20 


0.50 


0.53 


1.45 


-0.36 



Table 3.1: Coefficients A and a for modified depth calculation [17] in longitudinal 
development case, equation 3.1. 
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Figure 3.1: Longitudinal development of gamma ray induced showers with E t h 
lMeV particle detection threshold. 
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3.2 Lateral Development of Extensive Air Show- 



ers. 

As an extensive air shower cascades through the atmosphere ultrarelativistic parti- 
cles and photons are produced mainly in the forward direction. However, because 
electrons and positrons suffer multiple Coulomb scattering off the electric fields of 
nuclei and photons undergo Compton scattering off atomic electrons, the particles 
are spread out and the shower attains a lateral distribution. The density of parti- 
cles is greatest near the shower core, the trajectory that the primary gamma-ray 
would have had if it did not interact with the atmosphere. The average number 
of photons p(r, t, E ) per unit area at a distance r from the shower core agrees 
well with the Nishimura-Kamata-Greisen (NKG) formula, with modified depth 
parameter [17]: 

p(r,t,E )^^^-f(r/r m ,s 2 ) (3.2) 

where 

\ \ / r \ S2-2 / r \ S2-4.5 

f(r/r m ,s 2 ) = — -— — — 1 + — 

2tt 5(s 2 ,4.5 - s 2 ) \r m J \ r m J 



t 2 = t + b(E th ) 

s 2 = 

t 2 + 2y 

with B(x,y) being beta function so that 2ir J °° f(r/r m , s 2 )(r /r m )d(r /r m ) = 1, 
f'm — being Moliere scattering unit, r rn = = 9.7g ■ cm~ 2 or about 110 meters 

at the elevation of Milagro (E s = m e c 2 4rr / a « 21MeV) and N(E ,t) is total 
number of photons at the depth t [20] given by equation 3.1. 

It was also found in [17] that the lateral density distribution of electrons is 
well represented by the same expression if the scattering unit is substituted by 
r' m = r m /2. (This fact that shower photons are spread farther from the shower axis 
than the electrons is a consequence of the fact that photons do not lose energy by 
ionization and can travel larger distances than electrons.) The values of parameters 
b(E t h) are given in the table 3.2. 
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Eth, 


electrons 


photons 


MeV 


b 


b 


1 


0.45 


0.83 


5 


-1.22 


-1.49 


10 


-2.57 


-3.45 


20 


-4.22 


-5.51 



Table 3.2: Coefficients b for modified depth calculation [17] in lateral distribution 
case, equation 3.2. 

The average density of photons and electrons per unit area as a function of 
distance from the shower core is illustrated on figure 3.2 

3.3 Temporal Distribution of Extensive Air Shower 
Particles. 

Once an air shower develops lateral structure, one can speak about the shower front 
- the forward edge of the advancing cascade. The arrival time T of the earliest 
particle hitting a plane perpendicular to the shower axis at a distance r from the 
core provides the information concerning the shape of the front. Since the electron 
component of the lateral distribution is attained due to multiple scattering, it is 
lower energy electrons and positrons which propagate farther away from the axis. 
These travel at lower speeds and one expects them to be delayed with respect to 
the energetic one's at the core. The photons are expected to be more prompt than 
electrons, but are also delayed with respect to core due to greater distance traveled. 
The average arrival time of the first particle as a function of the core distance for 
both electron and photon components is illustrated on figure 3.3. 

It appears that the front of each of the components assumes a parabolic shape. 
The fluctuations of T around its average appears smaller for the photon component 
than for electron one [17] and, in general, are quite small. The thickness of the 
shower is defined by the arrival time distribution of particles with respect to the 
front at the given distance r, and is, therefore, due to lower energy electrons and 
photons originating later in the shower cascade. 
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(a) Lateral distribution of electron compo- (b) Lateral distribution of photon compo- 
nent, ncnt. 

Figure 3.2: Density of particles as a function of core distance for gamma ray 
induced showers with E th = lMeV particle detection threshold. Curves are nor- 
malized to the total number of particles in the respective components. Plots are 
made for showers at the depth of 20 radiation lengths. 
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Figure 3.3: Average arrival time of the first particle as a function of the core 
distance. 
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3.4 Cherenkov Radiation. 



One of the methods of detection of high energy electrons and positrons and there- 
fore of the air showers is with the help of Cherenkov radiation. Cherenkov radiation 
arises when a charged particle traverses a dielectric medium with a velocity v which 
is greater than the speed of light in the medium c/n. Here n is the index of refrac- 
tion of the medium which, in general, depends on the wavelength of the emitted 
radiation and c is the speed of light in vacuum. The radiation is emitted by atoms 
and molecules polarized by the moving particle. According to Huygens' principle, 
the partial waves will interfere to create the total wavefront propagating at an 
angle 9 = arccos — with respect to the velocity of the particle. At lower speeds, 
some energy of the particle is still lost due to polarization of the medium, however 
no Cherenkov radiation results. 

The condition for existence of radiation — = 4- < 1 can be expressed in terms 

nv pn 1 

of the particle's energy: 



mc 2 



E > = (3.3) 

V 1 ~ h 

where m is the rest mass of the particle. For example, the index of refraction of 
water is n water = 1.35, and therefore the threshold energy for electron to produce 
Cherenkov radiation (and therefore detection threshold) is equal to 0.759 MeV. 
For a muon it is 0.157 GeV and for proton 1.40 GeV. 

The number of Cherenkov photons produced per unit path length by such a 
particle with charge Ze and per unit wavelength [18]: 



d 2 N 2naZ 2 



1 - 



dxdX X 2 V P 2 n 2 (X)J 

or, in terms of emitted energy: 

d 2 E _ ixe 2 Z 2 ( _ 1 \ 
dxdX ~ e A 3 y 1 ~ (3 2 n 2 {X) ) 

2 

where a = 47T e £ohc is fine structure constant. Here, we have made explicit the de- 
pendence of the index of refraction on wavelength. Cherenkov radiation is emitted 
preferably in the short wave region (blue/violet end of the visible spectrum). 
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Note that for ultrarelativistic particles (3 ~ 1, the Cherenkov angle in water is 
approximately 42° and the emitted energy becomes independent of the particle's 
energy. This fact can be used to provide absolute energy calibration of photode- 
tectors. 



3.5 Detection of Extensive Air Showers. 

Air showers produced by high energy photons consist mainly of electrons, positrons 
and lower energy photons, and thus hint at methods of detection. If the energy of 
the primary photon is greater than about 10 TeV, then there are enough particles 
in the cascade reaching the surface of the Earth to enable the detection by arrays 
of scintillation counters where energy of charged particles is converted into flashes 
of visible light. These flashes are detected by photoelements. The detectors of 
this type are called extensive air shower arrays (EAS arrays). CASA-MIA and 
Tibet are examples of such arrays. The determination of the arrival direction and 
of the primary energy of air showers, sampled by a detector array, makes use of 
the lateral distribution function over a wide range of distances: the arrival time of 
the shower particles is used to determine the direction, the measurements of the 
particle density as a function of core distance together with the direction provide 
information about the energy of the primary gamma ray. 

Cascades, initiated by a primary photon with energy below several hundred 
GeV do not reach the ground, but can be detected by Cherenkov radiation that 
charged particles emit in the air as the cascade develops. Such detectors are 
termed air Cherenkov telescopes, such as Whipple and HEGRA. Detectors of this 
type, being optical devices, make use of the fact that most of the Cherenkov 
light is emitted in the forward direction of the primary particle and rely on their 
angular resolution — truly a telescope type of a measurement. Arrival direction 
and energy of the primary particle is inferred on the basis of the imaged longitudinal 
development of the air shower. 

Other types of detectors which were used to detect secondary particles are 
bubble chambers, spark and ionization chambers. Such methods are rarely used 
nowadays. Instead, atmospheric Cherenkov telescopes and extensive air shower 
particle detector arrays constitute the two major ground-based techniques. 
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3.6 Cosmic Rays. Difference Between Cosmic 
Ray and Gamma Ray Induced Showers. 



The discussion has been concerned so far with the air showers produced by pri- 
mary gamma rays. However, among the particles entering the Earth's atmosphere 
gamma rays present a fraction of less than 0.1%. About 79% of the particles are 
high energy protons and 14% are alpha particles. The rest are nuclei of heavier 
atoms: carbon, nitrogen, oxygen, iron,... 

Just as in the case of gamma rays, high energy hadrons also initiate air showers 
when they enter the atmosphere. However, the processes in a hadronic cascade are 
quite different from the electro-magnetic one. 

In such a cascade an incident hadron undergoes strong interactions with the air 
nuclei in which protons, neutrons, mesons and hyperons are created in quantities 
determined by their relative cross-sections. The most numerous particles are that 
of the pion triplet. 

The charged pions have a lifetime of about 2.6 • 1CT 8 seconds and, if they did 
not interact first, decay into muons and neutrinos: 

7T + -> /i + + 
7T~ -> IT + 

which gives rise to a muonic component of the cascade. If the energy of charged 
pions is high, then, due to relativistic time dilatation, they will have the oppor- 
tunity to interact with a nucleus rather than decay, thus producing secondary 
particles which replenish the hadronic component of the cascade. Neutral pions, 
on the other hand, have a much smaller lifetime of about 10 -16 seconds and decay, 
dominantly, into photons (it ^7 + 7). If energy of photons is above critical, they 
will initiate electromagnetic cascades. It is due to presence of the electromagnetic 
component the hadron induced air showers are very similar to those initiated by 
gamma rays. Muons, produced in the cascade, have a lifetime of 2.2 • 10~ 6 seconds 
and high energy ones survive to the sea level because of time dilatation. 

Thus, in a hadron cascade, the secondary particles either interact again, decay 
or are absorbed as a result of ionization energy loss. The cascade builds up to 
a shower maximum as in an electromagnetic cascade after which the numbers 
decrease. Because decay muons do not interact by strong interactions and loose 
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less energy in bremsstrahlung radiation then electrons, they primarily loose energy 
by ionization and therefore have a slower decrease. Because detection of high 
energy gamma rays relies on registration of secondary particles, the extensive air 
showers produced by cosmic rays constitute background noise for ground-based 
gamma-ray detectors. Special techniques and algorithms have to be developed to 
suppress this noise in order to increase the sensitivity to photon primaries. The 
presence of muons in the hadron induced shower is often exploited. 

3.7 Milagro: a Next Generation of EAS Array. 

The design of a detector has to reflect the goals and the obstacles mentioned above. 
The great success of air Cherenkov telescopes is due to their excellent angular 
resolution and good background rejection. With advantages come its limitation. 
Because the energy threshold for electrons to produce Cherenkov radiation in air 
is about 21 MeV (n a i r = 1.000293, formula 3.3), it cuts the number of detectable 
particles by a factor of 2 (table 3.1, equation 3.1). This bounds the energy of the 
primary gamma rays detectable by the technique to be greater than about 100 GeV. 
Furthermore, the telescopes are narrow-field-of-view optical devices, therefore, they 
can observe only a small portion of the sky at a time during cloudless, moonless 
nights. This severely impacts on their ability to detect transitory sources and 
perform observations of extended sources. 

Extensive air shower arrays, on the other hand, can observe the entire overhead 
sky and can operate 24 hours a day regardless of weather conditions. However, 
scintillation counters typically cover only a small fraction of area compared to 
lateral extent of the shower, and therefore detect only a small fraction of secondary 
particles leading to higher energy threshold, typically above 10 TeV. 

Milagro achieves an energy threshold of 400 GeV by using water as a detec- 
tion medium and by locating the detector at relatively high altitude. Using water 
allows Milagro to detect not only charged secondary particles of the air shower, 
but also secondary photons via Cherenkov radiation of cascades initiated by them 
in the water. This is important because secondary photons present a large frac- 
tion of particles reaching the ground level (table 3.1, figure 3.1). Also, the time 
distribution of the foremost particle appears to be narrower for the photon compo- 
nent than for the electron one [17] rendering the efficient conversion of secondary 
photons into Cherenkov radiation by the water medium even more essential. 
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Figure 3.4: Schematic diagram of the Milagro detector. 



Milagro is a covered, light tight pond filled with purified water and instrumented 
with photo detectors. Relativistic, charged particles produce Cherenkov radiation 
as they traverse the water. The radiation is emitted in a cone-like pattern with an 
opening angle of not greater than 42°. This allows Milagro to instrument a large 
surface area of the detector with a sparse array of photo detectors. Methods such 
as these have been used for decades in high energy physics [19]. 

As shown in figure 3.4, Milagro has been configured with two planes of photo 
detectors. The upper layer is used to measure the air shower front, providing the 
information needed to reconstruct the primary particle direction. The lower layer 
is used to detect penetrating muons or hadrons and aid in rejecting cosmic ray 
induced showers. 

Milagro fits the classification of extensive air shower array, shares many prin- 
ciples but improves on detection methodology. 
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Chapter 4 

The Milagro Detector 



4.1 Physical Components of the Detector. 

The Milagro detector is located in the Jemez mountains, about 45 km west of Los 
Alamos, New Mexico, at -106°40'37" East longitude, 35°52'43" North latitude. 
The central detector consists of 723 photomultiplier tubes (PMT) deployed in a 
60 m x 80 m x 8 m pond filled with water. The elevation of the reservoir above sea 
level is 2650 meters, which translates to an atmospheric overburden of 750 g/cm 2 . 

The buoyant photomultiplier tubes are held below the water surface by anchor 
cords, attached to a grid of weighted PVC pipes. The spacing of the grid is 
2.8 meters. The length of each string was calculated for each PMT so that the 
PMTs would all lie in a horizontal plane and form a two layer structure. Each 
PMT is floating upright with its photocathode facing upwards. Each PMT is also 
surrounded by a conical baffle to block internally reflected light. 

The 450 PMTs of the top layer are deployed under 1.4 meters of water and are 
used to measure the arrival direction of the shower. Two hundred seventy-three 
additional PMTs are located near the bottom of the pond under 6 meters of water 
and are used to distinguish photon- and hadron-induced air showers. The top layer 
is called the shower layer and the bottom one — the muon layer. 

The complete set of PMT string attachment points was surveyed after the grid 
was installed. Vertical positions of the PMTs were verified after the pond was filled 
with water. The final accuracy of PMT coordinates is estimated to be ±0.03 m 
in horizontal and ±0.01 m in vertical directions. Coordinates of the PMTs are 
used for shower reconstruction and detector calibration. Because PMTs are used 
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to detect Cherenkov light produced by air shower particles traversing the water, it 
was necessary to block the external light from entering the pond with an opaque 
cover and to provide for high transparency of the water. The 1 mm thick cover 
was made of two black layers of polypropylene with an internal polyester scrim. 
The cover can be inflated to allow access into the pond. High quality of the water 
was achieved with continuous recirculation and filtration of the water by a set of 
carbon filter, 1 — /im filter, ultraviolet lamp and 0.2 — /im filter. The attenuation 
lengths of the water in the recirculator and that of the pond are indistinguishable 
and equal to 13.4 ± 0.5 meters [22]. 

A lightning protection system was installed around the experiment as a safe- 
guard against possible strikes. 

4.2 The photomultiplier tube. 

The photomultiplier tube is a vacuum device used to transform very faint light 
signals into electric ones. It consists of a photocathode, a set of dynodes and an 
anode. A photon incident on the photocathode causes the emission of an electron 
(called photo- electron (PEJ) into the vacuum tube via the photoelectric effect. 
This electron is directed towards the first dynode by the focusing fields. When it 
hits the dynode, secondary electrons are emitted, each of which is guided towards 
the next dynode. The dynodes are kept at different electric potentials and thus 
cause acceleration of electrons to sufficient energies for secondary emission to take 
place. As this cascade develops, the number of electrons increases exponentially 
and by the time it reaches the anode the number of electrons is about 10 7 for every 
electron emitted at the photocathode. 

Besides amplification, the following points motivate the choice of the Hama- 
matsu R5912 SEL PMT model: 

• Cherenkov radiation generated by shower particles in the water is emitted 
mostly at the violet end of the visible spectrum. The selected model is 
sensitive in the 300 - 650 nm wavelength range. 

• Quantum efficiency, the ratio of the number of PE produced at the cathode 
to the number of incident photons, is rather high: 0.2-0.25 for the selected 
model. 
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• Time resolution. Time resolution of a PMT is limited by the fluctuations in 
the cascade development starting from the photocathode all the way down to 
the anode. The variance of this transit time increases as the input intensity 
decreases and is 2.7 ns at the lowest input (1 PE produced). 

• Late pulsing. Late pulses on the output of a PMT are believed to occur 
when electrons are being reflected off the first dynode and produce secondary 
electron only upon second re-entry into the dynode assembly caused by the 
focusing fields. The probability of late pulses decreases for high light inten- 
sities because the probability that all electrons produced at the cathode are 
reflected decreases as number of them increases. 

• Pre-pulsing. Pre-pulses are considered to be caused by the light hitting the 
first dynode directly, thereby producing an anode signal which precedes the 
main one induced by the photo-electron. The Probability of pre-pulsing is 
higher for larger light intensities. 

• After-pulsing. After pulses are thought to be caused by residual gas molecules 
in the tube being ionized by electrons and subsequently hitting the photo- 
cathode and producing electrons. The signal induced by these electrons 
follows the main anode one. 

• Saturation. Saturation is the effect of decrease of the amplification coefficient 
of the PMT for high light intensities. This is caused by the inability of the 
dynodes to accelerate the increased number of the secondary electrons to 
sufficiently high energy. 

Overall, the characteristics of the selected PMT model were found suitable for 
the application. 

4.3 Electronics. Time Over Threshold. 

Each PMT is connected to its own electronic channel with a single cable used to 
supply high voltage for the tube and to deliver a signal pulse from it. The signal 
carries two pieces of information about the Cherenkov light hit: the time when the 
hit occurred and it strength. 
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Figure 4.1: Time over threshold. 
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The time of the hit is determined by the time when the amplified signal crosses 
the discriminator with the preset threshold. This time is termed start time. A 
time-to-digital converter (TDC) is used to digitize the start time and to pass it to 
the data acquisition computer. Except for the saturation, the output anode charge 
of the PMT is proportional to the number of Cherenkov photons which struck the 
tube, therefore, the output charge is the measure of the hit strength. 

Conventionally, analog-to-digital converters (ADC) are utilized to digitize the 
charge information. Because of cost, speed and dynamic range limitations of typical 
ADC's, a time-over-threshold (TOT) method was adopted instead. The idea of 
the method is simple (figure 4.1). A capacitor C is quickly charged by the PMT 
signal with total charge Q and then is slowly discharged via a resistor R. The 
voltage on the capacitor as a function of time is given by: 

V(t) = %e-*l RC 

If the time during which the voltage exceeds the preset threshold is measured, 
then it is seen that 

time over threshold ~ In Q 

Thus, time over threshold is related to the initial charge Q and provides a way 
to measure the signal strength. Implementation of this method allows the use 
of a single multi-hit TDC module to digitize both start time and the time over 
threshold counterparts. 

In practice, two discriminators with different threshold levels (low and high) 
are used to guard against pre- or after-pulsing of the PMTs. These accompanying 
pulses are usually small and do not cross the high threshold. Therefore, the use 
of high threshold information is preferred where available. The low threshold was 
set to about 1/4 of the average signal produced by a single photoelectron and high 
threshold to about 7 PE. 

Each shower layer PMT signal crossing the low discriminator also participated 
in the trigger formation process by generating a 25 mV, 300 ns long pulse [21]. 
The analog sum of all such pulses is sent to a discriminator with a threshold 
corresponding to 60 PMT being hit. If the discriminator is fired the trigger pulse is 
generated. Setting the trigger threshold much higher than 60 PMT would increase 
the energy threshold of the detector. Making it much lower would enable single 
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Figure 4.2: Schematic diagram of a PMT channel. 



muons to trigger the detector. The choice of 300 ns window is motivated by the 
duration of near horizontal shower propagation through the detector. The trigger 
provides a common stop for all TDC's, all start times are referenced to it. 

The absolute time of each triggered event is provided by a Global Positioning 
System clock. A schematic diagram of a PMT channel is presented in figure 4.2. 



4.4 Timing Edges. 

Time-over-threshold pulses generated by both discriminators (if crossed) are mul- 
tiplexed into a series of time edges. (Figure 4.3) It is the polarity 1 and time of 
these edges which is being recorded as raw data from a PMT. Strong hits cross 
both thresholds and result in a 4-edge event, weak ones cross low threshold only 
and produce 2 timing edges. 

In practice many edge events can be observed. One edge events can result 
when the other edge (or edges) is truncated by the trigger time window. Four edge 

1 Polarity is the direction in which the threshold was crossed: going up or down. 
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Figure 4.3: Timing edges. 



events can be a sequence of two 2-edge events. Multi edge hits are also possible. 
A filtering algorithm had been designed to clean the hits. It bases its decisions on 
the polarity of hits and their relative time separation. On average, about 8% of 
PMT hits are rejected as having invalid edge sequences. The accepted PMT hits 
are characterized by start time (LoStart) and time over low threshold (LoTOT) 
and their high threshold counterparts (HiStart, HiTOT) if available. 



4.5 Monte Carlo Simulation of the Detector. 

The simulation of the detector response is done in two steps: (1) initial inter- 
action of the primary particle (proton or photon) with the atmosphere and the 
subsequent generation of secondary particles, and (2) detector response to the 
secondary particles reaching the detector level. The CORSIKA [23] air-shower 
simulation code provides a sophisticated simulation of secondary particle develop- 
ment in the Earth's atmosphere induced by a primary particle with energy up to 
10 20 eV. Within CORSIKA, the VENUS code is used to treat hadron-nucleus and 
nucleus-nucleus collisions at high energies whereas GHEISHA is used at low en- 
ergies (< 80 GeV). Electromagnetic interactions are simulated using EGS 4 code. 
The atmosphere adopted in CORSIKA consists of nitrogen, oxigen and argon in 
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volume proportions of 78.1%, 21.0% and 0.9% respectively. The density variation 
of the atmosphere with altitude is modeled by 5 layers. In CORSIKA a flat atmo- 
sphere is adopted which is a good approximation up to zenith angles of about 70° 
where discrepancy with the shperical one reaches one radiation length. 

The simulation of the detector itself is based on GEANT [24]. All of the sec- 
ondary particles in a shower cascade reaching the Milagro are used as input to 
the GEANT with uniform random placement of the shower core around the detec- 
tor. GEANT simulates the electromagentic and hadronic interactions of particles 
in the pond and results in the set of simulated times and pulse heights at each 
PMT. This information is saved in the same format as the real data and is used 
to establish properties of the detector. The detailed simulation of the electronics 
and phototubes continues to be the area of active reseach. Further improvements 
and understanding are expected. 
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Chapter 5 
Calibration. 



5.1 System Setup and Goals. 

The calibration system has been designed to reflect the physics goals of the detector 
and is used to obtain parameters needed to transform the raw counts to physically 
meaningful arrival times and light intensities which then can be used for event 
reconstruction. Despite the considerable effort that has been made to construct all 
PMT channels of the detector as uniformly as possible, in order to achieve the high 
precision required for the event reconstruction the remaining variations between 
channels have to be compensated for. A separate set of calibration parameters is 
determined for each PMT channel. 

The desire to determine the positions of events on the Celestial Sphere with 
systematic errors much less than the expected angular resolution (which is about 
1°) dictates that the locations of the photo-tubes be known to about 10 cm accuracy 
in horizontal direction and 3 cm in vertical, and PMT timing accuracy to about 
1 ns. Photographic and theodolite surveys were used to ensure accurate PMT 
position determinations. 

To achieve the stated time accuracy it is important to calibrate the TDC con- 
version factors, compensate for pulse amplitude dependence of TDC measurements 
(known as slewing correction) and synchronize all TDCs (find TDC time off-sets) 
to the required accuracy. Time over threshold to photo-electron conversion must 
be determined to convert all PMT amplitude measurements to a common unit 
for each event. All of the above is achieved with the help of the laser calibration 
system. 
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Figure 5.1: Calibration system setup. 



The calibration system is based on the laser — fiber-optic-diffusing ball concept 
used in other water Cherenkov detectors [25]. A computer operated motion con- 
troller drives a neutral density filter wheel to attenuate a 300 picosecond pulsed 
nitrogen dye laser beam. The selected dye emitted light at 500 nm. The beam is 
directed to one of the thirty diffusing laser balls through the fiber-optic switch (see 
figure 5.1). Part of the laser beam is sent to a photo-diode. When triggered by the 
photo-diode, the pulse-delay generator sends a trigger pulse to the data acquisi- 
tion system. A laser fire command is issued by the motion controller, providing full 
automation of the calibration process. The balls are floating in the pond so that 
almost every PMT can be illuminated by more than one light source. Complete 
description, operation, analysis, problems and suggestions are described in [29]. 
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5.2 Timing Calibration. 



5.2.1 Slewing Calibration. 

Time slewing is the dependence of the reaction time of a PMT together with 
its electronics on the intensity of the incident light. The amount of light in the 
calibration system is regulated by the filter wheel whose transmission property 
ranges from completely opaque to completely transparent thus enabling the study 
of the effect. The ideal situation of no slewing should reveal itself as independence 
of PMT registration times with respect to photo-diode when intensity of light in 
the calibration system is changed. 1 

In practice, it takes longer for weak pulses to cross a discriminator threshold, 
and thus results in a delay. This is illustrated on the figure 5.2. 

Because the light intensity is characterized by time over threshold, the amount 
of slewing is studied with respect to this variable. The method of generation of 
such a slewing curve is explained with the help of the figure 5.3. 

This diagram is high/low threshold independent, the scheme is applied inde- 
pendently for both threshold levels. The slewing curve is a plot of T start vs ToT. 
Note, that T start includes the propagation time of pulses in the water and in the 
optical fiber, delays associated with the details of the particular PMT channel and 
common detector trigger delays. Because common offsets are irrelevant, after cor- 
recting for relative fiber delays and water propagation times which effectively shifts 
the curve up or down, this will become a complete timing calibration. The slewing 
correcting curve is found by fitting a polynomial to T start vs ToT. If a PMT gets 
a slewing curve from more than one laser ball, the curve resulting from the max- 
imal light illumination is chosen. An example of a slewing curve is presented in 
figure 5.4. The TDC's used in Milagro employ a common stop, thus larger values 
of T start correspond to earlier times. 

5.2.2 TDC Conversion Verification. 

The time to digital converters measure time in the units of "counts" . According to 
LeCroy 1887 FASTBUS TDC specifications one count corresponds to 0.5 nanosec- 
onds. This was verified by insertion of known variable delay in the photo-diode 

1 The photo-diode is thought as having no slewing problem because intensity of light incident 
on it is constant. 
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Figure 5.2: Slewing and Time over Threshold. The weaker the pulse (or, equiva- 
lently, the higher the discriminator threshold) the latter the registration time is. 
Infinitely large pulse will cross the threshold without any delay, or zero slewing. 
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Figure 5.3: Conceptual drawing of a measurement performed for timing calibration. 



36 



100 200 300 400 500 600 700 800 900 1 000 



100 200 300 400 300 600 700 800 900 1000 



(a) Example of the T sta rt vs ToT data 
points. 



(b) Example of the polynomial fit. The 
values in each bin are the means of the 
Tstart obtained by a fit of a Gaussian 
distribution in the corresponding ToT 
bin. 



Figure 5.4: Plots showing a typical data obtained from the calibration system 
the purpose of timing calibration. Units of both axis are TDC counts. 
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trigger chain with the help of DG535 digital delay/pulse generator by Stanford 
Research Systems. The shift in the start times (T start ) of all PMTs was found to 
be very consistent with the 2 count per 1 ns conversion. This assured that all TDC 
clock speeds are the same and meaningful interpretation of time across all PMT 
channels is available. 

5.2.3 Fiber-Optic Delays. Speed of Light in Water. 

The last step of timing calibration is correcting the slewing curve for both the light 
propagation time in water between the PMT and the laser ball used to generate the 
curve and the relative delay in the fiber attached to that ball. In order to correct 
for the propagation time, coordinates of PMTs and balls and the speed of light 
in water have to be given. PMT and laser ball coordinates are known from the 
survey. (If laser ball coordinates are deemed to be inaccurate, they can be found 
from the calibration data: see [26, 27].) Because only a typical index of refraction 
of water is found in reference tables and because fiber-optic delays vary from ball 
to ball, these parameters have to be evaluated on the basis of the calibration data 
itself. 

The method to solve the problems uses the ability to cross calibrate a PMT 
using several laser balls. If Ti and T 2 are slewing corrected times registered by a 
PMT from two laser balls (1 and 2), then the difference (7\ — T 2 ) does not depend 
on delays in the PMT electronics (common to both measurements), but instead 
reflects the difference in propagation times from the corresponding balls and the 
relative fiber delays. Consider the following quantity: 

T T\ T2 + Ajjfrer -|- ^propagation 

If speed of light and fiber delays are correct, then r is zero within errors of 
measurement. When for given laser ball pair many PMTs with their r's are con- 
sidered, one sees that for PMTs located approximately half way between the balls 
^■propagation ~ and deviation of average r from zero is due to fiber delay difference. 
On contrary, PMTs located in the close proximity of either ball provide informa- 
tion about the speed of light: the usage of incorrect speed of light will widen the 
distribution of r's. Thus the analysis of all r's for all PMTs and laser ball pairs 
enables the determination of speed of light in water and relative fiber delays. 

Slewing curves are shifted by the appropriate propagation time and fiber delay 
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to yield the final timing calibration for each PMT. 

5.3 Photo-Electron Calibration. 

5.3.1 Occupancy Method. 

Time over threshold to photo-electron calibration is based on the well-known oc- 
cupancy method, which was one of the methods of PE calibration of Milagrito [28] 
and other water Cherenkov detectors [25] . The data used for PE calibration is the 
same as for the timing calibration which is obtained with the laser light passing 
through the filter wheel with different transparency settings (which are changed 
by rotating the wheel). The main task of the PE calibration is determination of 
the number of photo-electrons for a given incident light intensity (ToT). The basic 
method consists of two steps: 1. at low light levels calibrate using the occupancy 
method, and 2. extend to higher light levels using given attenuation properties of 
the filter wheel. 

1. Calibration at low light levels. 

At low light levels it is possible to measure the number of photo-electrons as 
function of ToT directly. The method is founded on the following assump- 
tions: 

• The number of photo-electrons produced at PMT's photocathode obeys 
a Poisson distribution: 

w = 

where A is the mean number of the photo-electrons produced. (This is 
justified because the probability of emitting a single photo-electron does 
not depend on the fact of possible emission of other photo-electrons.) 

• A constant light intensity source is used in calibration. 

Then, the probability rj that at least one photo-electron was produced while 
the photocathode was illuminated (the probability that the PMT "saw" the 
illumination) is called occupancy and is given by: 



39 



rj = p(n > 0) = 1 - P x (n = 0) = 1 - e~ 



Occupancy can be easily measured if the PMT is illuminated several times 
with the same pulse intensity: 

number of observed pulses 

V = 

total number of pulses 

Therefore, if n is known, then A = — ln(l — n). This method is applicable at 
relatively low light levels (when A < 2) because at high levels the occupancy 
saturates to 1 and a small measurement error in n will lead to a big error in 
A: 

AA = —^—An = e x An 
1 — n 

2. Calibration at high light levels. 

Calibration at high light levels is achieved by noting that in the absence of 
PMT/electronics saturation the number of photo-electrons is proportional 
to the light intensity at the photocathode. If the transmittance of the filter 
wheel T is known, then: 

A = aT 

where a is some coefficient which is different for each laser ball PMT pair. 
This coefficient can be found at low light levels where A is also known. Error 
of this method will grow linearly with light intensity, not exponentially as 
in low light level case. Thus, given transmittance T of the filter wheel, 
ToT-to-PE conversion can be achieved at any light level. When, in reality, 
saturation is present prescribing the number of photo-electrons using this as 
a requirement allows protection against the effect. Traditional amplitude- 
to-digital conversion methods which make use of the air shower data are 
susceptible to such non-linearity. 
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5.3.2 In Situ Filter Wheel Calibration. 



Traditionally, filter wheel calibration is obtained either from the manufacturer or 
by a separate measurement. In Milagro, in situ filter calibration method was 
developed which uses the same calibration data. The idea is grounded on the 
supposition that for any two sufficiently close levels of transmittance of the filter 
wheel (Ti and T2) there exists a PMT in the pond for which the occupancy method 
is valid on both light intensities. If Ai and A2 are corresponding photo-electron 
measurements, then: 

T2 _ A2 
T\~ AT 

This line of arguments is used to relate T 3 to T 2 , T 4 to T 3 and so on, leading 
to restoration of the transmittance levels for the entire wheel. Because absolute 
transmittance of the wheel is irrelevant one can always set T — 1. 



5.3.3 Dynamic Noise Suppression. 

Small amounts of radioactive elements present in the water, ambient light, Che- 
renkov light produced by the shower particles or thermo-electron emission cause 
signals on the output of PMT channels which constitute noise hits. The presence of 
these hits will lead to overestimation of occupancy implied by the laser light level 
and thus will damage the accuracy of the calibration. Dynamic noise suppression 
is a technique allowing the solution of this problem on the tube by tube basis. 

An event, registered by a PMT could have come from a laser pulse or from 
a noise hit which are independent, therefore the probability of observing a hit 
P(any) is given by: 

P[any) = P[laser + noise) = P(laser) + P(noise) — P[laser) • P[noise) 

where P {laser) is probability of observing a pulse due to laser (true occupancy), 
P [noise] is probability of observing a pulse due to noise. P[any) is measured by 
sending laser pulses into the pond, P[noise) is obtained by sending "random" (no 
laser light) triggers to the data acquisition system, then the occupancy is given by: 

, P[any) — P[noise) 
n = PUaser) = — i — / \ — - 
' V ; \-P[noise) 
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Figure 5.5: Filter wheel calibration curves obtained with different laser ball - 
PMT combinations are presented. Bank 2,3 represent laser balls 11-20 and 21- 
30 respectively. AS, MU — PMT's of air shower and muon layers respectively. 
Calibrations with ball 19 and muon layer PMT's are also shown with and without 
the noise suppression to illustrate its importance. 

This is indeed how occupancy was determined for the filter wheel calibration 
and for the ToT-to-PE conversion. To demonstrate the importance of the suppres- 
sion refer to figure 5.5 which shows filter wheel calibration curves obtained with 
and without noise suppression. 



5.4 Slewing Extrapolation. 

The maximum light level at which a PMT was calibrated depends on the distances 
to the laser balls, laser light output and the calibration system optics alignment. 
The range of ToT within which a tube is calibrated may be too narrow to allow 
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interpretation of the strongest hits generated by the shower particles. The desire to 
make use of these hits requires extrapolation of the calibration curves beyond the 
measured points by some plausible functional form. The functional shape depends 
on the PMT discriminator threshold level, amplification coefficients and so on. 
Instead of trying to make physical model of the channel a statistical approach is 
followed. 

It is believed that all PMT channels were manufactured to meet common char- 
acteristics. Therefore, the study of the channels' responses (calibration) can be 
viewed as a multiple (about 700 times, broken PMTs do not count) measurement 
of a single function of ToT. The fact that the calibration curves obtained for 
different channels are slightly different can be attributed to the "manufacturing 
imperfections" (such as unavoidable uncontrollable spread of characteristics of elec- 
tronic components) and can be treated as such. Thus, the curves for all PMTs can 
be viewed as particular realizations of the same random function [30] . 

Inasmuch as a random function is a mathematical concept of great complexity 
and in most general case it can be regarded as a non-denumerable set of scalar 
random variables, it is natural to seek the expression of the random function in 
terms of simpler random concepts: ordinary scalar random variables. Therefore, 
the following representation of the calibration curve X(t) is tried (here t or t' 
denote ToT): 

n 

X(t) = m(t) + £ V v <f> v 

Here, {<p V) v — 1, ..,n} are some non-random functions, {V u } are random vari- 
ables. This representation is said to be canonical if m(t) is the mean function 
(m(t) = M[X(t)]) and V^s are uncorrelated with zero mean [31]: 

M[K] = 0, D[V V ] = M[(K) 2 ] = D v 

M[K-V;] = n±v 

When such a representation is found, each particular calibration curve x(t) can 
be extrapolated using the expansion and the set of coefficients v v can be calculated 
in the range where calibration data is available. The construction of the canonical 
representation makes use of the average function m{t) and the correlation function 
K(t,t') = M[X°(t) ■ X°(t')}, where X°(t) = X(t) - m(t). It can be shown that 
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satisfy the conditions of the canonical representation if coefficients are cho- 
sen such that 



which is always possible if the number m of sample points {t k } is large enough 
(m > (n — l)/2) [32]. The root mean square deviation of the representation from 
the realization is given by: 



This scheme is used for extrapolation of slewing curves of timing calibration. It 
is implemented to the first order of the canonical expansion [32]. The only sample 
point t\ is the highest value of ToT available from the calibration data. Therefore, 
the extrapolation is obtained by: 



The typical extrapolation range needed to interpret shower data is less than 
100 ns in ToT leading to the error of the extrapolation being less than 0.7 ns [32]. 
Comparison of the extrapolated curves and measured ones obtained in different 
calibration runs yielded the measured extrapolation error of 0.55 ns [33] in good 
agreement with the expectations. 
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Chapter 6 

Event Reconstruction. 



6.1 Primary Particle Direction Reconstruction. 

The data registered in each event is used to infer the characteristics of the primary 
particle that created the shower: its direction of incidence, energy and type. 1 In 
order to be able to make the inference about the properties of the primary particle 
based on the observed quantities one has to assume the relation between the two. 
One of these relations is expressed by equation 3.2, the lateral distribution of the 
shower. If the arrival direction is known, then the depth t of the atmosphere tra- 
versed by the shower can be determined from geometrical consideration and then, 
equation 3.1 together with equation 3.2 relate the measured lateral distribution 
p{r) with energy of the primary particle and shower core position. From the prac- 
tical point of view, the fluctuations in the shower development, small size of the 
detector and fluctuations in its response make this program infeasible [34]. The 
detector is currently being upgraded with the array of water tanks which will aid 
in core location, allowing achievement of 30% energy resolution [34]. 

The task of determination of the primary particle direction is much simpler and 
is of highest importance (from the gamma ray astronomy point of view). The basic 
idea is as follows. Let us assume that the shower front is a plane, perpendicular to 
the direction of motion of the primary particle. Then, as in any air shower array, 
each detector % in the set records the local arrival time Tj of the shower front (see 
figure 6.1) which has to be contrasted with the model: a plane defined by a normal 

1 The PMT registration times and corresponding light intensities deduced from the raw counts 
and calibration parameters present the event data. 
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Figure 6.1: Schematic diagram of a shower array detectors and how timing infor- 
mation can be used to determine the direction of the primary particle ft. Lines 
represent measured arrival times, filled circles — the locations of detectors. 



vector n = (a, b, c), the direction of the primary particle: 

axi + by { + czi 
f v a z + b z + c z 

where, (xi,yi,Zi) are coordinates of PMT i, vq — is the speed of the shower 
front in air, T — is a common time offset for the whole event. Noting that Milagro 
is a flat array, thus Z; L = 0, Vi, and introducing notations: 



A = 



t>oV 'a 2 + b 2 + c 2 
we obtain the final model 



B 



vqV o? + b 2 + c 2 



Ti = Axi + Byi + T 

When coefficients A and B are found, the direction of the incoming shower is 
reconstructed as 



cos(0) 
sin(0) 
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VA' 2 +B 2 

9 = arcsin^ov 7 ^ 2 + B 2 ) 
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where (0, 0) are zenith and azimuth angles, and v , the speed of the shower 
front in air, is approximately equal to the speed of light. 

The method of determination of coefficients A and B from the measured times 
Ti is a x 2 fit where weights are chosen based on the recorded light intensity of the 
PMTs. This is because more energetic particles appear to suffer less fluctuations 
due to the combined effect of the shower development and PMT time resolution 
and thus are given higher weight. The weights were derived by studying the dis- 
tributions of Tj's [35]. These distributions are not Gaussian, as assumed by the x 2 
fit, therefore, several iterations are made. The experimental points are also occa- 
sionally just way off (such as times of incidental hits). Points like this are called 
outliers. They can easily turn a x 2 fit on otherwise adequate data into nonsense. 
Their probability of occurrence in the assumed Gaussian model is so small that the 
estimator is willing to distort the whole curve to try to bring them, mistakenly, 
into line. Therefore, on the first iteration only PMTs with hits of greater than 
2 PE are used. On the next iterations the PE restriction is gradually relaxed so 
that PMTs with weaker signals are allowed to participate. At the same time, the 
definition of outliers becomes more stringent. After the first iteration points with 
distance of more than two standard deviations away from the plane are removed 
from the fit on the second pass. Before the last, fifth, pass outliers are defined as 
all hits with distance of more than 0.5 standard deviations away from the plane. 
The results of the last iteration are recorded as the reconstructed direction of the 
primary particle. The number of PMTs that participated in the last pass are also 
noted as Nf it . 

More sophisticated algorithms were also attempted, they, however, required 
a significantly higher computer power without improving angular resolution by a 
noticeable amount [36] 

6.2 Core Location. 

The inability to make precise determination of the shower core impacts not only 
energy estimation but also the direction reconstruction. Indeed, as was mentioned 
in section 3.3, the shower front is not flat, but rather presents a paraboloid surface. 
Therefore, in order to be able to use the described plane fit algorithm, the shower 
front has to be "unfolded" into the plane. The amount of this curvature correction 
is certainly core distance dependent. Failure to apply curvature correction will 
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result in a systematic error in the direction reconstruction which increases with 
core displacement from the pond. Because there are many showers with different 
core positions, this will lead to degradation of average angular resolution. 

Thus, for the purposes of angle determination, the following core locating algo- 
rithm was adopted [37]. First the direction to the core from the center of the pond 
is determined by calculating the location of the provisional core. The provisional 
core position is found in two iterations. Iteration one: the average PMT coordi- 
nate is calculated where the weight is taken to be y/PE of the corresponding PMT. 
The second iteration is the same as the first, but where only PMTs within ±45° 
sector of the direction found on the iteration one are used. After the provisional 
core is found, the decision is made whether the shower core is on the pond or it 
is outside of it. The decision is made on the basis of profile of the average PE 
versus distance from the center of the pond for PMTs within ±60° sector of the 
provisional core direction. If it is declining sufficiently fast towards the edge of the 
pond, the core is adopted to be on the pond. Otherwise, it is off. On-pond cores 
are found as PE weighted average of PMT coordinates in the circle of 8 meters 
around provisional core. Off-pond cores are placed at a distance of 50 meters from 
the center of the pond in the direction defined by the provisional core. Only PMTs 
of the shower layer are used in this algorithm. The value of 50 meters is chosen 
because simulations indicate that this is the most probable core displacement for 
showers triggering the detector. 

In the early part of the detector running, a more primitive version of the algo- 
rithm was used. 



6.3 Curvature Correction. 

The origin of the curvature correction is traced back to the desire of using the 
plane model of the shower front in the direction reconstruction and to the details 
of the shower development underlined in section 3.3: the shower front is parabolic. 
It was shown that electron and photon components of the shower have different 
shapes of the corresponding fronts (figure 3.3). This means that the shape of 
the front, as sampled by the array, depends on energy threshold of the detectors 
and on the relative efficiency for photon conversion. For this reasons, methods of 
determination of the curvature correction from the real data were developed [35]. 
A particular representation of the correction was chosen and then parameters were 
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optimized in a multistage iteration process. The correction was shown to improve 
angular resolution compared to a plane fit. 

6.4 Sampling Correction. 

The curvature correction and core finding algorithms provide for primary particle 
direction reconstruction given that the times of the arrival of the shower front are 
measured by the detectors. However, as was mentioned in section 4.2, each PMT 
has quantum efficiency of about 0.2 — 0.25. This means that the first particle, 
carrying information about the shower front, is detected with some probability p. 
If this particle is missed, the tube is presented a chance of detecting the second 
particle in the shower and so on up to the total number of particles p(r)dxdy which 
make up the thickness of the front at the specific core distance r. Thus, if f n (T, r) 
is the time distribution of the n-th particle, f±(T,r) being the shower front, then 
the distribution g(T,r) of the detected arrival times is given by: 

g(T, r) = pf ± (T, r) + (1 - p)pf 2 (T, r) + . . . + (1 - p) n ~ l pf n {T, r) + . . . 

If nothing else is known, the g(T,r) has to be interpreted as the shower front 
instead of f\{T,r). This would certainly lead to poorer angular resolution. 2 The 
Milagro analysis tries to exploit the correlation that particles trailing the front 
are of monotonously decreasing energy, making it possible, effectively, to estimate 
the order number n of the detected particle. With regard to the core distance 
dependence, each of the f n (T,r) is narrower where density of particles p(r) is 
higher, that is towards the core. Thus, the sampling correction is aimed at inferring 
the arrival time of the front needed for direction reconstruction, it is due to finite 
detection efficiency of the PMTs and it is a function of core distance and light 
intensity registered by the PMT. Such estimated front arrival times are subject to 
higher fluctuations compared to direct measurements and therefore are given less 
weight in the direction reconstruction algorithm [35]. 

2 Particles following one after another have their arrival times ordered: T\ < < . . . < T n < 
. . .. Among these times, T\ obeys distribution fi, T2 obeys and so on where corresponding 
variances are believed to increase with n. Therefore the variance of g(T,r) satisfies inequality 
var{g{T 1 r)) > var(fi(T,r)) with equal sign if and only if g(T,r) = fi(T, r) — 100 % detection 
efficiency. 
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6.5 Cosmic Ray Rejection. 



The identification of the type of the primary particle (photon or cosmic ray) in 
Milagro makes use of the presence of muons and hadrons in cosmic-ray initiated 
cascades. In the top layer of the detector their presence is obscured by the general 
illumination by the shower particles, in the bottom layer, however, they may be 
noticed. Penetration of one of these particles to the muon layer should lead to 
clusters of high light intensities. In the case of a primary photon, the illumination 
of the bottom layer should be uniform. It was found with the help of Monte Carlo 
simulations that the parameter X 2 defined as 

PEmax 

is sensitive to the type of the primary particle. Here N 2 is the number of 
PMTs in the muon layer with registered light intensities greater or equal to 2 
PE, PEmax — is the maximum intensity detected in the event by the muon layer 
([38] and [39]). The simulated distributions of parameter X 2 for gamma- and 
proton- induced showers are presented on figure 6.2. The X 2 distribution for the 
proton dominated data (also shown on figure 6.2) agrees with the simulated proton 
distribution well. This match is a check of the Milagro Monte Carlo simulations. 

It is clear that events with large values of the parameter are more likely to be 
gamma ray induced. It was determined that the optimal cut for gamma rays is 
accepting events with X 2 > 2.5. At this cut, about 90% of protons are rejected 
while about 50% of gamma rays are accepted. 

Attempts to develop a more complicated criterion have not yet yielded a sig- 
nificant improvement over the present scheme. 
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1 2 3 4 5 6 7 8 9 v 10 

Figure 6.2: Distributions of parameter Xi- Dashed line is Monte Carlo simulations 
for photon induced showers, dotted line is simulation for proton induced showers, 
solid line is the data. Each distribution is normalized to 1. 
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Chapter 7 

Performance of the Milagro 
Detector. 

7.1 Operation of the Milagro. 

The Milagro detector operates largely unattended in a reliable and stable manner. 
Automatic alerts are generated under serious error conditions such as loss of elec- 
trical power, an abnormal event rate, or overheating of the electronics. The nature 
of the error, time of the day and weather conditions determine the response time 
and restart of the data taking. Less serious problems can be corrected remotely. 
There are also scheduled down times to accommodate repair activities and cali- 
brations. Other than that, data are acquired continuously. During operation of 
the detector some PMT-electronics channels cease to work and have to be turned 
off. This problem has typically been traced to water leakage of the under-water 
connectors linking coaxial electric cable with the PMT. Repairs of these connectors 
performed once a year are a major part of the scheduled maintenance down time. 
The repair time is chosen to be September when the pond water is the warmest to 
facilitate the scuba diving to retrieve PMTs. 

7.2 Angular Resolution. 

The cosmic ray shadow of the Moon has been used to measure the angular reso- 
lution of an EAS array above 50 TeV [40]. At TeV energies, the geomagnetic field 
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will displace the Moon cosmic ray shadow. Another estimate of the angular resolu- 
tion is obtained by studying A eo — the difference between directions reconstructed 
by the same "color" PMTs if the detector is imagined to be painted in the white 
and black squares of a checkerboard. (The PMT numbering scheme was designed 
in such a way that the color of the PMT is defined by the parity (even/odd) of its 
number.) The quantity A eo is not sensitive to certain systematic effects, common 
to the both parts of the detector, however, in the absence of these effects the dis- 
persion of A eo is twice the overall angular resolution [41]. From the studies of the 
Moon shadow, the A eo and the Monte Carlo simulations, the estimated angular 
resolution of the instrument including pointing effects is 0.75° for all events with 
N fit > 20. 

7.3 Absolute Energy Scale. 

The absolute energy scale can be determined by examining the magnetic displace- 
ment of the shadow of the Moon. A preliminary estimate of the median energy of 
the detected events is 640±70 GeV which is in excellent agreement with simulations 
which predict median energy of 690 GeV. 

7.4 Cosmic Ray Rejection at Work. 

All observations to date indicate that the flux of gamma rays from the Crab nebula 
is constant which makes it a very useful source for testing the sensitivity of different 
instruments. The cosmic ray rejection method was tested on the Crab. If the data 
is analyzed without application of the Xi cut, the Crab is observed with significance 
of 1.4. Analysis of the data passing the rejection criterion yields significance of 
5.4. (The notion of significance is discussed in section 8.2.) The data set used in 
this study covers the period between June 8, 1999 and April 1, 2002 [42]. 
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Chapter 8 

Data Analysis Technique. 

8.1 Coordinate Systems on the Celestial Sphere. 
8.1.1 Equatorial Coordinate System. 

Inasmuch as distances to the majority of astronomical objects are much bigger than 
the size of the Earth and its orbit, it is sufficient to specify directions to the objects 
regarding distances as equal and infinitely large. Within this framework, the stars 
appear to be located on the surface of an imaginary sphere with the observer at its 
center. This sphere is called the Celestial sphere. Various coordinate systems on 
the sphere are used, all defined by the corresponding choices of reference circles. 
Horizon and equatorial coordinate systems are important examples employed in 
this work [43]. Figure 8.1 illustrates the following definitions. 

C - the observer. 

CP - line parallel to the axis of rotation of the Earth. 
P - North celestial pole. 

Celestial Equator - is defined as intersection of a plane perpendicular to CP at 
point C and the celestial sphere. (This is the projection of Earth's equator.) 

Z - the zenith, intersection of the celestial sphere with the outward continuation 
of the plumb line at the observer's location. 

IPCZ = 90° — (f)- definition of <fi. Angle <fi is astronomical latitude of the observer 
on the Earth. 
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Figure 8.1: Definitions of horizon and equatorial coordinate systems. 
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Horizon - is defined as intersection of a plane perpendicular to CZ at point C 
and the celestial sphere. 

N,S - are North and South of the Horizon. N and S are defined as the intersection 
of a great circle PZ centered at observer C with the Horizon. Arc PZLS is 
called local reference celestial meridian. 

X - is a star. 

LZCX =ZX= z - zenith distance of star X. 

A - azimuth, is a dihedral angle between reference meridian and ZXC plane. 

LLCX =LX= H - hour angle, is a dihedral angle between reference celestial 
meridian and that of the star PXD. 

IX CD =XD= 5 - is Declination of a star (Dec). 

The equatorial coordinate system of hour angle and declination (H, S) is built 
around the axis of rotation of the Earth whereas the horizon one of zenith and 
azimuth (A, z) uses a plumb line as the reference. The law of cosines for trihedral 
angles applied to the spherical triangle XPZ (figure 8.1) two times yields the 
relation between the systems: 

cos(90° — 5) = cos z cos(90° — 0) + sin z sin(90° — 0) cos A 
cosz = cos(90° - 5) cos(90° - 0) + sin(90° - 5) sin(90° - 0) cos# 

or, 

sin S = sin cos z + cos sin z cos A 

sin z sin A 
tan if = ■ ■ - 

cos cos z — sin z sin cos A 

The value of observer's latitude defined as the complement of the angle LPCZ 
is assumed to be known. Both of these coordinate systems are local coordinate 
systems in the sense that both of them revolve with the Earth. If the object X 
is stationary in space, due to the Earth's rotation it will appear to be moving 
and its local coordinates will be changing. With this motion, the declination 5 
of the stationary object in the equatorial system will remain constant, while its 
hour angle H will change incrementing by 360° when the Earth makes one full 
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revolution around its axis CP. The time required to complete such a revolution 
is called sidereal day. In contrast, the solar day or universal day is defined as the 
time between two appearances of the Sun on the local reference celestial meridian. 
The universal time is different from the sidereal time due to Earth's orbital motion 
around the Sun. 

Equatorial celestial coordinate system is defined as declination 5 and right as- 
cension a of the object, where a = H r — H, H r is the hour angle of the vernal 
equinox. This coordinate system is independent of the Earth's rotation and ob- 
server's longitude. Local sidereal time is defined as H r expressed in the units of 
time. 

Local sidereal time as well as the geodesic latitude of Milagro are provided by 
the GPS system clock. In general, there is no exact relation between astronomical 
latitude defined above and the geodesic one, but according to indications in litera- 
ture [44] the difference is of the order of 15" in magnitude, much smaller than the 
angular resolution of the detector and thus can be neglected. 

8.1.2 Galactic Coordinate system. 

The reference plane of the galactic coordinate system is the disc of our Galaxy 
(i.e. the Milky Way) and the intersection of this plane with the celestial sphere 
is known as the galactic equator, which is inclined by about 63° to the celestial 
equator. Galactic latitude, b, is analogous to declination, but measures distance 
north or south of the galactic equator, attaining +90° at the north galactic pole 
(NGP) and -90° at the south galactic pole (SGP). The galactic latitude of the 
star X on the figure is arc YX and is north (figure 8.2). 

Galactic longitude, /, is analogous to right ascension and is measured along the 
galactic equator in the same direction as right ascension. 1 The zero-point of galac- 
tic longitude is in the direction of the Galactic Center (GC), in the constellation 
of Sagittarius; it is defined precisely by taking the galactic longitude of the north 

1 This sense of rotation, however, is opposite to the sense of rotation of our Galaxy. The Sun, 
together with the whole Solar System, is orbiting the Galactic Center at the distance 28,000 light 
years, on a nearly circular orbit, moving at about 250 km/sec. It takes about 220 million years 
to complete one orbit (so the Solar System has orbited the Galactic Center about 20 to 21 times 
since its formation about 4.6 billion years ago). Considering the sense of rotation, the Galaxy, 
at the Sun's position, is rotating toward the direction of I = 90°, 6 = 0°. Therefore, the galactic 
north pole, defined by the galactic coordinate system, coincides with the rotational south pole of 
our Galaxy, and vice versa. 
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Figure 8.2: Definition of Galactic coordinate system in its relation to the equatorial 
coordinate system. 

celestial pole to be exactly 123°. The galactic longitude of the star X on the figure 
is given by the angle between GC and Y. 

The galactic north pole is at RA = 12:51.4, Dec = +27:07 (2000.0), the galactic 
center at RA = 17:45.6, Dec = -28:56 (2000.0). The inclination of the galactic 
equator to Earth's equator is thus 62.9°. The intersection, or node line of the two 
equators is at RA = 18:51.4, Dec = 0:00 (2000.0), and at I = 33°, 6 = 0° [45]. 
Transformation between celestial coordinate system and galactic one are given, 
therefore, by [46]: 

cos b cos(/ - 33°) = cos 5 cos(a - 282.85°) 

cos6sin(Z - 33°) = cos S sin(a - 282.85°) cos 62.9° + sin 5 sin 62.9° 
sin b = sin 5 cos 62.9° - cos 5 sin(a - 282.85°) sin 62.9° 

8.2 Significance of a Measurement. 

According to its definition [47], a statistical hypothesis is a theoretical prediction 
or assertion about the distribution of one or several measurable quantities. It is 
therefore desired to devise a rule of checking whether the data are consistent with 



58 





(a) View of the Galaxy from the North 
Galactic pole. Galactic longitude grid is 
shown. The Galaxy, at the Earth's po- 
sition, is rotating toward the direction of 
I = 90°. 



(b) Side view of the Galaxy. Galactic lat- 
itude grid is shown. The arrow indicates 
the direction of increment of Galactic lon- 
gitude which is opposite to the Galactic 
rotation. 



Figure 8.3: Artist view of the Milky Way Galaxy from aside. 
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the hypothesis, called null hypothesis H . Such a rule, however, can not be hoped 
to tell whether the hypothesis is true or false, it can only insure that in the long run 
of experience wrong propositions will not be accepted too often. A naive method 
of testing is to calculate the probability that a certain character x of the observed 
data would arise if Hq were true. If the probability is small, this is considered 
as an indication that the hypothesis is false and it is rejected. However, if x is a 
continuous variable, then probability of any value of x is zero. Therefore, tests are 
based on the notion of critical region. The critical region w is the subset of space 
W of all possible values of x such that x falling into w leads to the rejection of 
the hypothesis. The hypothesis is accepted or remains in doubt in all other cases. 
It follows then, that if there are two alternative tests for the same hypothesis, the 
difference between them consists in the difference in critical regions. It also follows 
that H can be rejected when, in fact, it is true {error of the first kind)] or it can 
be accepted when some other alternative hypothesis H\ is true {error of the second 
kind). (Admittance of the existence of an alternative hypothesis is clear, otherwise 
the null hypothesis would not be questioned.) The probabilities of errors of these 
two kinds depend on the choice of the critical region. The critical region is said 
to be the best for testing hypothesis H with regard to Hi if it is the one which 
minimizes the probability of the error of the second kind among all regions which 
give the same fixed value of the probability of the error of the first kind. The 
construction of the best critical region, resulting in the most efficient test of H 
with regard to Hi was considered by Neyman and Pearson in [47]. The problem 
is solved for the general case of simple hypotheses. A hypothesis is said to be 
simple if it completely specifies the probability of the event; it is composite if the 
probability is given only up to some unspecified parameters. 

Critical regions w(e) corresponding to different probabilities e of errors of the 
first kind are engineered before the test is performed. When experimental data x 
is obtained, the smallest e is found such that x G w(e). It is then said, that based 
on the experimental data, the null hypothesis can be rejected with significance e. 



A typical astronomical experiment consists of two independent observations of 
some regions of the sky yielding two measured counts Ni and N 2 made during 
time periods ti and t 2 respectively with all other conditions being equal. The null 
hypothesis being tested is that Ni and N 2 constitute a sample of size 2 from a 
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single Poisson distribution (adjusted for the duration of observations) which is due 
to common illumination. The alternative is that the two measurements are due to 
Poisson distributions with unrelated parameters. In other words, the alternative 
hypothesis Hi is of the form: 

while the null hypothesis is also of the composite type, obtained from the Hi 
by setting Ai = A2: 

that both measurements are due to the same unknown count rate A = Ai = A2. 

Despite the fact that po(Ni,N 2 ) admits the existence of regions similar to the 
sample space W with regard to the parameter 2 A, no common best critical region 
with the respect to every alternative admissible (Ai 7^ A2) hypothesis pi(Ni,N 2 ) 
can be found. 

Therefore, the following practical method is adopted which is to consider 

jVi — a No 

S = j= =Sf= , a = W*2 (8.1) 
V Ni + a z N 2 

Arguments can be made (appendix A) that for large N\ and N 2 , under the null 
hypothesis, S obeys the standard normal distribution in the vicinity of zero. The 



2 According to publication [47], the region w is called region of size e similar to sample space 
W with regard to parameter (3 if 



/ 



Po(x, (3)dx = e = const 



for all values of parameter [3. The hypothesis po(Ni,N 2 ) satisfies the special case considered 
in [47], for which 

9 4= A+B * 

where cj> = —7^, coefficients A and B may depend on f3, but not on x. For the case of interest, 



po dp 

= X, ( j > =^±^-2,A = -l,B- T 
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bounds of this approximation — So < S < S enlarge when smallest of the TV's 
increases (equation A.l). The critical regions w(e) with respect to the alternative 
hypothesis Ai > A 2 are, therefore, defined as S > S(e), where: 

1 r 00 , 
e = . — / e 2 ax 

V27T JS(e) 

and provide the way to calculate the significance of the measurement e. (The 
critical regions w(e) with respect to the alternative hypothesis Ai < A 2 are defined 

as S < S(e), e = -j= jf^ e~^dx.) Because of the one-to-one correspondence 
between e and S(e), the significance of a measurement can be quoted in the units 
of S (see table 8.1). The error on the value of e due to this approximation is less 

than -l-f£e-4dx. 



\S(e)\ 


e 


1.0 


0.15866 


2.0 


0.02275 


2.5 


6.21 • 10~ 3 


3.0 


1.35- 10" 3 


3.5 


2.326 • lO" 4 


4.0 


3.17- 10~ 5 



Table 8.1: The correspondence between l^l and e for the alternative hypothesis 
Ai > A 2 or Ai < A 2 [48]. 



8.3 Effective Area. 

Considering the detector and the atmosphere together as parts of a giant appa- 
ratus, it is obvious that not all particles entering the atmosphere will trigger the 
detector. (Particles on the other side of the Earth will not be registered by the 
detector.) Thus, it is meaningful to speak about detection probability, or more 
precisely, about conditional probability for a particle to be detected given its type, 
its direction of incidence, its energy and displacement of its trajectory from the 
detector's center (core position). Then, the number of detected events is given by: 
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N obs = Y,J P(detect\k, E, 0, x, y) F{k, E, 0) dE dQ dx dy 

k 

where F(k, E, 0) is the signal function, the total number of particles of k-th 
kind emitted with energies between E and E + dE in the directions -j- © + dQ 
of local sky and landing in the area dxdy during the measurement period. Note, 
that the signal function does not depend on the core position (x,y), which is a 
reasonable assumption: landing coordinates on the Earth have nothing to do with 
the emission process in the sky. Therefore, the integration over core positions can 
be performed and the result is called the effective area: 

/+oo 
P(detect\k,E,Q,x,y)dxdy = A eff (k,E,Q) (8.2) 
-oo 

N obs = W A eff (k, E, Q)F(k, E, 0) dE d& (8.3) 
k J 

The term effective area is used because A e ff(k,E,Q) has units of area and 
can be interpreted as the geometrical area of a fictitious detector with detection 
efficiency of 100% independent of a particle's core position over its extent and zero 
outside. Because air showers may have lateral extent of 300 meters or more, they 
may trigger the detector even if the core is very far from the detector. That is why 
effective area can be much larger than the geometrical area of a real detector. Due 
to the fact that the signal function F(k, E, 0) and the number of observed events 
N b s are related via equation 8.3, effective area becomes an important characteristic 
of a detector. Detectors of relatively small sizes can be calibrated using test beams 
of particles of interest and the effective areas can be measured experimentally. 
For ground based detectors such as Milagro where the Earth's atmosphere is a 
component of the detector the effective area has to be computed with the help 
of computer simulations. Definition 8.2 of effective area provides a direct method 
for its determination: Monte Carlo integration. To do this, N to tai showers are 
simulated with core distances chosen uniformly over a sufficiently large area A 
and those which trigger the detector are counted. 

Ae " {k ' E ' 0) - N total (k,E,e) Ao 
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Such Monte Carlo integrations were performed for photon and proton primaries 
with energies in the range of O.lTeV-lOOTeV and zenith angles between 0° and 
45°. Appropriate tables are saved and used in the analysis. (Dependence of the 
effective area on the azimuth angle is disregarded in this calculation, therefore, in 
what follows the effective area is denoted as A ljCr (z, E).) It is important to note 
that because the atmosphere can be considered as one of the constituents of the 
apparatus, changes of its conditions can lead to variations of the effective area. 



8.4 Determination of Gamma Ray Flux from a 
Source. 



When the presence of a source is established by means of the hypothesis test of 
section 8.2, it is desirable to measure the features of the source, the parameters of 
its signal function. 

In the case of a point source, the signal function has the form 



where <&(E,t) is the differential flux of the emitted particles 
([&(E, t)] = m 2.s ec . eV ) and z(t) is the trajectory of the source in the local sky 
(z is the zenith coordinate). If the source is steady, <£> does not depend on time 
and the signal function becomes: 



where T(z) is the effective amount of time spent by the source in the zenith 
region between z and z + dz. The case of an extended source may be reduced to a 
collection of point sources, and if it is uniform, it will be absorbed by T(z). The 
physical processes responsible for the gamma ray emission from the Galactic plane, 
as well as those of the cosmic rays are believed to be steady and are believed to 
produce power law differential flux type spectrum: 




observation period 




F(E,z) = $(E)-T(z) 
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ME) = |^**(> E ) ■ E-°» 



E 



leading to the signal function prototype: 

F k (E, z) = ^^F k (> Eo) ■ E~^ ■ T(z) (8.4) 

where F k (> E ) is the integral flux of particles with energies greater than E , 
F k (> Eq) = § k (E)dE. Then, the number of gamma ray and cosmic ray induced 
events from the source region can be calculated with the help of equation 8.3: 

N, = ^^tF,(> E ) [ [ A,(z,E) E~ a ~< T(z) dz dE 



N cr = ^ c _ r Qcr+ 1 1 Fcr(> E ) J J A cr (z,E) E~ a ~ T(z) dz dE 



Then the ratio of the integral gamma and cosmic ray fluxes from the source is 



related to the ratio -^ 2L , the experimentally available quantity by: 



F 7 (> E ) 1 K 



7 



For 

(> E Q ) r?(a 7 , a cr ) N c , 



(8.5) 



where 



< \ a -'- 1 1 1 M*> E ) E ~ a ~< T{z) dz dE 

r?(a 7 ,a cr ) = E J 



a cr - 1 u JJ A cr (z, E) E~ a ° r T(z) dz dE 

The dimensionless coefficient 7?(a 7 , a cr ) is the energy and transit averaged ratio 
of effective areas for the given source. The effective areas were obtained with the 
Monte Carlo simulation, however, the ratio of them 77 is believed to be less sensitive 
to possible imperfections of the simulations than effective areas themselves. 

Formula (8.5) is the principal formula for interpretation of measurements for 
the assumed signal function (8.4). 
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Chapter 9 



Background Estimation. Null 
Hypothesis. 

9.1 Time Swapping Method. 

Most air showers detected are produced by charged cosmic rays that form a back- 
ground (chapter 3) to the search for gamma initiated showers from a source. Be- 
cause of their charge and because of the presence of random magnetic fields in the 
interstellar medium, the cosmic ray particles lose all memory of their initial direc- 
tions and sites of production, and can be regarded as forming isotropic radiation. 
Emission from a gamma ray source would appear as an excess number of events 
coming from the direction of the source. Therefore, a search for emission is a sta- 
tistical test with the null hypothesis that there is no source conducted by means 
of two observations. Indeed, one of the measurements could be direct counting of 
the events in the angular bin in celestial coordinates containing the source and the 
other observation, called off-source, providing information about background level 
in the neighborhood of the source bin. If according to the results of the test the 
two measurements are inconsistent with each other, it is said that the presence 
of the gamma ray emission is established. Otherwise, the measurements do not 
contradict to the null hypothesis that the observations are due to isotropic back- 
ground and no source detection can be claimed. This, however, does not preclude 
its existence. In other words, chapter 8 is the navigation map for source detection: 
identify the candidate source by its coordinates on the sky, perform measurements 
and calculate significance, if significant, estimate of the source function parameters 
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can be tried. 

Special consideration must be given to the significance test (section 8.2) because 
of the ambiguous statement that "two independent observations ... with all other 
conditions being equal" . The off-source observation can be performed at the same 
time as the on-source one utilizing the wide field of view of the detector, or it can 
be performed at a different time making measurement in the same directions of the 
field of view. (Due to the Earth's rotation, the off-source bin may present itself in 
the directions of local coordinates, previously pointed at by the source bin.) Both 
of these stipulations contradict to the conditions of "being equal": if observations 
are done at the same time, then non-uniformity in the acceptance of the array to 
air showers due to detector geometry must be compensated for; if observations are 
done at different times, then uniform operation of the detector must be assured. 
The mechanism of such an equalization is called background estimation. 

The widely accepted method of background estimation [49, 50] follows the 
second path by recognizing that no major changes in the detector configuration 
are usually made on the short time scale and by taking advantage of the rotation 
of the Earth. Therefore, and if the cosmic ray background is due to isotropic 
radiation, the number of detected events as a function of local coordinates x and 
time t can be written in the form: 



Here R(t) is overall event rate, G(x) - - acceptance of the array such that 
5 field of vie W G(, x ) dx = 1. The number of background events in the source bin, is 
then given by 



where 4>(x, t) is equal to zero if x and t are such that it translates into inside 
of the source bin, and is one otherwise. 

In the context of section 8.2, the number of events N s obtained from the ob- 
servation of the source celestial bin can be considered as Aq and number N b of 
estimated background events can be considered as aN 2 . Then equation 8.1 can be 
written as: 



dN(x,t) 



G{x) ■ R{t) dxdt 



(9.1) 





S = 



N s -N b 



(9.3) 
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where a is the ratio of effective exposures of the two measurement as before. 

The task now is to evaluate integral in (9.2). Recall [52] that Monte Carlo inte- 
gration of f v g(x)f(x)dx can be performed by generating x uniformly distributed 
over the area of integration V and then calculating the sum: 

r V N 

/ g(x)f(x)dx= —22g(xi)f(xi) 

jv iv i=1 

where N is the sample size generated. However, if f(x) can be interpreted as 
probability density function, then integration can be performed by generating x 
according to distribution f(x) and calculating sum: 

1 N 

g(x)f(x)dx = -T7^2g(xi) 

JV iv i=1 

The proof is evident if the integration variable is changed to y — ff^ f(x)dx, 
y G [0, 1]. Therefore, if N is total number of events detected during integration 
time of equation 9.2, N — J R(t)dt, then, introducing r(t) = R(t)/N , both G(x) 
and r(t) can be interpreted as probability density functions in corresponding spaces 
and integration of (9.2) can be done by means of Monte Carlo: 

AT N 

^ = TFE(1-<K^)) (9-4) 
iV i=i 

where (x,t) are distributed according to joint probability density G(x)r(t). A 
list of all coordinates of the detected events is regarded as a sample from the G(x) 
distribution, while list of all times — as the one from r(t). Therefore, sample from 
G(x)r(t) distribution can be generated from the data by randomly associating an 
event's local coordinate x with an event's time t among the pool of detected events. 
The so created coordinate-time pair is called a generated event. The accuracy of 
Monte Carlo integration increases as the square root of the number of generated 
events. 

In practice, Monte Carlo integration is performed by substituting each real 
event's arrival time by a new time from the list of registered times of collected 
events in a finite time window. That is why the method is referred to as time 
swapping method. The swapping is repeated (5 times per each real event, (3 typically 
being around 10. Then, Nj, is estimated as number of generated events in the source 
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bin divided by (3 (equation 9.4). The method entails that the second, off-source, 
celestial bin is defined by the regions of the sky which have the opportunity to 
present themselves into the same set of local coordinates as the source bin does 
due to the Earth's rotation during the time period of integration. 

9.2 Accounting for Signal Events During Swap- 
ping. 

Despite the fact that the time swapping method possesses the desired equalization 
properties, namely, the generated events all obey correct local angle distribution 
and correct timing distribution, even if the detector stability assumption (equation 
9.1) holds, the application of the method introduces difficulties. Indeed, equation 
(9.3) for significance of a measurement was obtained with the assumption of inde- 
pendence of Ni and N 2 - The time swapping realization of background estimation 
as described includes on-source events N s in the calculation of iV&. The problem 
can be brought to an extreme as in the case of sighting of the Polar Star. In this 
case, the source does not present any apparent motion in local coordinates be- 
cause it lies on the axis of rotation of the Earth. The off-source bin does not exist 
and measurement can not be performed. In the framework of the time swapping 
method, the background N b will be estimated to be exactly equal to N s report- 
ing zero significance. This is clearly unsatisfactory. There are two directions to 
proceed: one — modify the significance calculation to reflect the dependence of 
measurements, two — modify the time swapping method to regain independence 
while keeping all other properties intact. The latter one is followed in the current 
analysis. 

To regain independence of the two measurements, the events from the source bin 
should not participate in the background estimation. However, simply removing 
these events from the swapping procedure will destroy its foundation that the list 
of local coordinates and times represent samples from G(x) and R(t) respectively. 
This problem is solved by the following algorithm [51]. 

Denote by N out (x,t) the number of detected events originated outside of the 
source bin, and R ou t(t) their total event rate, then it is readily seen that 

dN 0Ut (x,t) = (f)(x,t)G(x)R(t) dxdt 
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Integrating this equation with respect to t and x the system of equations on 
unknown G(x) and R(t) is obtained: 




G(x) J <f>(x, t') ■ R(t') dt' 
R(t)f<i)(x',t)-G(x') dx' 



The numerical solution of these integral equations (see appendix B) provides 
R{t) and G(x) based on data N out (x) and R ou t{t) from the outside of the source 
bin. One can arrive at these equations by invoking the maximum likelihood prin- 
ciple based on the Poisson distribution of detected events with average fj,(x, t) = 
(f)(x,t)G(x)R(t)dtdx. 

The event rate R(t) is considered to be constant on the very short time scale 
(24 sidereal seconds amount to rotation of the Earth by 0.1°, much smaller than 
angular resolution of the detector which defines "very short") and therefore is 
saved as a histogram. This is the distribution of detected times estimated using 
off-source events only. Generated event times are drawn from it. The task now 
is to generate a sample from G(x) using off-source events. This sample should 
contain N(x) = G(x) J R(t)dt events with given local coordinates x, however, the 
number of off-source events available is N out (x) = G(x) J (fi(x,t)R(t)dt. Therefore, 
missing events are created by swapping available ones 



number of times. This can be done by choosing actual number of swaps from 
a Poisson distribution with parameter (1 + a'(x)). The described method is the 
Monte Carlo generation scheme applied to the equation (9.4) where only events 
from the off-source region are used. 

Two remarks are in order. First, the value a'(x) depends on local coordinate 
x. Second, a' has a similar meaning as a in equations (9.3) and (8.1): ratio of 
on and off source exposures. However, if the region of interest is a subregion of 
the excluded bin, then a in equation (9.3) pertains to the subregion, whereas, a' 
corresponds to the excluded bin. 

The importance of this modification is demonstrated on figure 9.1 where results 
of Monte Carlo simulations are presented. The figure shows excess number of events 
(7V 7 = N s — N b ) as a function of Galactic latitude before (figure 9.1(a)) and after 
(figure 9.1(b)) the change. The 25% signal loss is recovered by the modification. 



1 + a'(x) 



G(x)JR(t)dt 



G(x)JR(t)dt 
Nout(x) 



G(x) J (p(x,t)R(t)dt 
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(a) Excess number of events iV 7 as a func- 
tion of Galactic Latitude. Source region is 
not excluded. 



(b) Excess number of events iV 7 as a func- 
tion of Galactic Latitude. The region of 
±7° around Galactic equator is excluded. 
Modified time swapping method is used. 



Figure 9.1: Plots showing the results of Monte Carlo simulations with uniform 
Galactic signal flux being 0.0088 that of background in the region of ±5° around 
the Galactic equator. The expected bin content is about 205000. 
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9.3 Significance of a Measurement and the Time 
Swapping Method. 

Significance calculation according to equation (9.3) relies on the integral evaluation 
(9.2) of the estimated number of background events N b . The time swapping method 
employs Monte Carlo integration and thus introduces additional fluctuations in 
the estimate of iV&. The significance calculation has to reflect this fact and has to 
reduce to the equation (9.3) in the limit of infinite number of generated events, 
(3 ^ oo. Also, as was mentioned earlier, it is the time swapping method which 
defines the off-source region, therefore, the ratio of on- and off-source exposures a 
has to be calculated within the framework of the method and can not be supplied 
externally. In what follows, the explicit dependence on local coordinates x will be 
omitted and restored at the end. Also, new definitions are introduced, old ones 
are as before. 

Let N be the number of events available for swapping in the off-source region. 
iV is obtained as a result of a measurement and therefore is a sample of size one 
from a Poisson distribution with some parameter /i: 

P,(N) = f£ e -. 

These events participate in the swapping procedure, during which the total 
number of generated events M is drawn from a Poisson distribution with parameter 
(1 + a')[3N: 

P (M)- ((l + «OW M r -(iW)^ 

ni+a')my M > e 

During this Monte Carlo integration, m events out of M end up in the source 
region of interest which obeys a binomial distribution with some parameter p: 

B p (m,M)=CZp m (l- P ) M - m 
Combining all of the above, the probability distribution of m is given by: 

oo oo 

P(m) = ]T B p (m,M) ]T P {l+a , )pN {M) ■ P,(N) 

M=m N=0 
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Calculation of dispersion D{m) of m makes use of the following facts: 

D(m) = E[m 2 ] - E[m] 2 



N~P^(N): E[N}= fi, E[N 2 ]= fi + fi 2 

m ~ B p (m, M) : E[m\ = pM, E[m 2 ] = pM(l - p) + p 2 M 2 
and results 

E[m] = (1 + a')f3pn 
E[m 2 ] = (1 + a')Ppn + (1 + a') 2 /3V^(l + A*) 

Therefore, 

D{m) = (1 + a')/9pA* (1 + (1 + 

The unknown probability p(x) is estimated in such a way that the average E[m] 
is equal to the obtained: m — (1 + a')j3pfi. The parameter /i is estimated to be 
equal to the number of events available for swapping N out {x). Thus, the dispersion 
of the number of generated events in the source region is estimated as: 

s / m{x) \ 

Dim) = mix) 1 + 



by: 



N out (x)J 

Therefore, the dispersion of the number of background events Nb = ^m is given 



The ratio „ b ^\ is the measurement of a, the second term vanishes as 3 in- 

N out (x) ' i 

creases. Thus, finally, the significance of a measurement within framework of time 
swapping is given by 
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Six) = N a (x)-N b (x) = JU± 

y/N a (x)+a(x)N b (x) + ^N b (xy N out (x) 



9.4 Compounding Independent Experiments. 

When two independent observations of a source are made, a question of reporting 
total significance arises. The problem is addressed by noting that the two observa- 
tions imply existence of two sets of numbers (N^\ N b ^\ a^) and (N^ 2 \ N b 2 \ a^) 
needed for significance calculation in each of the experiments according to equation 
(9.3). 1 The combined significance is obtained by 

s {NU - NP) + (N?) - N®) 

COmP ° Und y/(NP + aWN b {l) ) + (nP + a^N b (2) ) 

This expression is a direct extension of equation (9.3) and bears all the prop- 
erties of significance. Examples of such a compounding are combination of inde- 
pendent experiments conducted at different times, or combination of independent 
but concurrent sightings of different fragments of a source. 

However, the two experiments can be combined in a different way, if the stability 
assumption (9.1) is common for both observations, namely, the two data sets can 
be united and the time swapping procedure can be applied to this common data 
set. Such a compounding may lead to increased total significance if otherwise 
useless information obtained during experiment 1 can prove to be valuable for 
experiment 2 and vise versa. The situation is illustrated on the figure 9.2 where 
two contiguous experiments are presented. The shaded area is the united off- 
source region, whereas, the heavily shaded is the combination of the independently 
useful off-source regions of the two experiments. It is seen that in united data set 

N (union) = N (l) + N {2)^ N (union) _ ^(1) + ^(2) ( statisticaUy equal), but a {uni ° n) < 

qX 1 . 2 ) because a's are ratios of exposures of on- and off-source regions, the latter 
one being expanded by the unification. Therefore, significance 



1 The discussion is of the general character, therefore, correction due to time swapping (equa- 
tion 9.5) is ignored in this section. 
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Field of view Hour Angle 



Figure 9.2: Illustration of difference between union and compounding of experi- 
ments. 
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jvr (union) _ at (union) 

q S b 



(union) ~ a ( union ) jy(union) 

is statistically higher than S compoun d because of its lower denominator. The 
increased sensitivity provides an incentive to broaden the limits of integration 
(9.2) for as long as stability assumption (9.1) allows. 



9.5 Detector Stability Assumption Test. Diurnal 
Modulations. 

Despite the fact that no reconfigurations to the detector on the short time scale 
are made, the acceptance of the array G(x) depends on transmission properties of 
the atmosphere and therefore, verification of the stability assumption (9.1) is still 
warranted [53]. 

The test of stability would be a comparison of two acceptances Gi(x) and 
G 2 {x) measured at different times t\ and t 2 . On physical grounds, the detector 
possesses a certain degree of azimuthal symmetry, so does the atmosphere, there- 
fore acceptance is considered as function of zenith and azimuth angles G(z, A). A 
measurement is a generation of histograms G\,{z,A) and G 2 (z, A) from the data 
for a certain duration of time around ti and t 2 . (This interval was chosen to be 30 
solar minutes.) These histograms are then to be compared using a significance test. 
It has to be recognized that presence of sources on the sky will mimic instability 
(section 9.7), therefore, zenith and azimuth angle distributions alone are compared 
instead. The test is based on ideas of compounding probabilities from independent 
significance tests [54] and is implemented as a series of x 2 tests of Gi(x) and Gj(x) 
(yielding x 2 (ti,tj)) an d then combined Xtotai(^) f° r time separation At = ti — tj 
is calculated: 

Xtotal(^)= J2 X 2 (U,t,) 
At=ti-tj 

The test statistic Xtotai (At) so obtained follows a x 2 distribution with m to tai = 
J2At=ti-tj m (ti,tj) degrees of freedom if observed differences are of random nature 
only. The average of Xtotai i s equal to m tota i while its variance is equal to 2m to tai- 
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Figure 9.3: Results of stability test on data collected between July 19, 2000 and 
September 17, 2000 are presented. Horizontal axis is time separation At in 30 
minute units, vertical axis is corresponding ^7 '/(At) • pl°ts are the test of 

zenith distributions, middle and bottom ones are the test of azimuth distributions 
for 0° < z < 30° and 30° < z < 50° respectively. On the left all data were used, on 
the right, only events passing Xi > 2.5 cut. Solid horizontal line is the expected 
value of one if the stability assumption holds. 
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The results of such a test are presented on figure 9.3 as the plot of total , , 4 . The 

1 to ^ m tota i(At) 

closer the plotted value to the expected one, the better the assumption of stability 
(9.1) holds. It is seen from the plots that the degree of violation of the assumption 
grows with time separation At as might be expected, but then it drops before grow- 
ing again. This is interpreted as presence of a periodic component which insured 
that two acceptances G\{z) and G2(z) separated by 24 solar hours are "closer" to 
each other than, say, those separated by only 12. Azimuthal distributions show vi- 
olation of the stability assumption to a lesser degree. This statement is considered 
to be true after application of the hadron rejection cut as well. Thus, despite the 
fact that no human intervention on the short time scale is made, the acceptance 
of the detector changes. As was seen in section 9.4, the duration of the validity of 
the assumption (9.1) defines the time integration limits in equation (9.2) which, in 
turn, defines the off-source region. Inasmuch as search for very weak signal requires 
sensitivity pushed to its limits, care is taken to improve the stability assumption 
in particular with regard to zenith angle dependence. 

The investigation of changes in zenith distribution (simplified by the above 
noted periodicity) was performed in the following way [55]. Every half hour, a 
zenith distribution is generated from registered events. These distributions are 
compared to the average distribution accumulated over a one day period. Half 
hour distributions and the average one are normalized and then subtracted to 
give a plot of distribution's shape change. It was observed that the shape of this 
change is approximately constant with amplitude varying from half hour to half 
hour. Therefore, the improved stability assumption is chosen to be of the form: 

N(x, t) = G{x)R{t)e e(t)K{z)+q ^ m) dt (9.6) 

where 6{t) is the amplitude of the correction at time t, K{z) is the polynomial 
correction function coefficients of which are obtained from the above study (see 
appendix C), q(0{t)) is the normalization factor so that / G(x)e e ^ K( -^ +q(e ^ dx = 1. 

The choice (9.6) of the form of time dependence of a distribution is motivated 
by the Darmois theorem [56]. According to this theorem, exponential class (9.6) 
of probability density functions is the only class which admits number of sufficient 
statistics for unknown parameter 6 independent of the number of observations. In 
particular, if (Zi, Z 2 , Zn) is a random sample of size N from the distribution 
(9.6), then 
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N 

Y = Y,K(Z i ) 



is sufficient statistics for 9. The importance of sufficiency is explained by the 
fact that Y exhausts all the information about 9 that is contained in the sample. 
By differentiating 



twice with respect to 9 we obtain the following two relations for expectation 
and variance: 



Both of these properties allow estimating parameter 9 from the data. Indeed, 
because 9 = corresponds q = (no correction to average distribution), both 9 
and q can be determined as 



Expectation and variance and therefore qi and qii are determined from the sam- 
ple and thus are known. The example of the average daily amplitude dependence 
is shown on figure 9.4. The value of the amplitude is typically within ±4 ■ fCT 5 
range. The plot can also be used to justify the choice of half hour intervals for the 
amplitude measurement. (This, however, is already seen from the figure 9.3.) 

9.6 Accounting for Modulations in the Frame- 
work of Time Swapping Method. 

Incorporation of the improved stability assumption (9.6) into the framework of time 
swapping is equivalent to modification of Monte Carlo integration where generated 




E[K(Z)} 



qr(6) 



var[K(Z)\ 



qff(9) 
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Figure 9.4: Average daily dependence of the zenith correction amplitude 9 derived 
from July 19 - September 17, 2000 Milagro data after application of the hadron 
rejection cut. 
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(a) Distribution of x 2 for 90 bin histogram. (b) Distribution of \ 2 f° r 360 bin his- 
Solid line is the best fit of \ 2 distribution togram. Solid line is the best fit of \ 2 
function with 84.11 degrees of freedom. distribution function with 337.0 degrees of 

freedom. 

Figure 9.5: Plots showing Monte Carlo simulated distributions of x 2 of the differ- 
ence of histograms with 90 and 360 bins. Both histograms have the same number 
of entries (5713). 

events are distributed according to G(x)R(t)e e ^ K ^ +q ^ d ^ instead of G(x)R(t). 
This is accomplished by means of the rejection method [52]. In this method, an 
event is drawn from the G(x)R(t) distribution as before. But now, the candidate 
event time is accepted only if a random variable p distributed uniformly between 
and P(z) = max {t} e e ® K W +q W» is such that p < e mK(z)+ q (e(t)) _ (p aramete rs 
9 and q{9) are assumed to be known at this stage.) If this inequality is violated, 
another candidate time is sought and then tried. This is continued until time is 
accepted, resulting in generation of event. Inasmuch as generation of the events is 
based on the improved stability assumption, the rest of the data analysis algorithm 
is preserved. 

In order to illustrate how well the correction works real and generated distri- 
butions are compared using a x 2 test. The zenith angle distribution is split into 
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(a) Distribution of \ 2 for standard sta- 
bility assumption (equation 9.1). Solid 
line is the best fit of \ 2 distribution func- 
tion with 93.84 degrees of freedom. The 
goodness of this fit is characterized by 
x 2 /ndf = 506.2/50. 



(b) Distribution of \ 2 f° r improved sta- 
bility assumption (equation 9.6). Solid 
line is the best fit of x 2 distribution func- 
tion with 84.10 degrees of freedom. The 
goodness of this fit is characterized by 
X 2 /ndf = 91.34/86. 



Figure 9.6: Plots showing distributions of \ 2 of the difference of simulated and real 
zenith angle distributions for every half hour with and without correction. The 
expected number of degrees of freedom is 84.25. Both histograms have the same 
number of entries (2073). 
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(a) Distribution of \ 2 f° r standard sta- 
bility assumption (equation 9.1). Solid 
line is the best fit of x 2 distribution func- 
tion with 339.6 degrees of freedom. The 
goodness of this fit is characterized by 
X 2 /ndf = 69.82/83. 



(b) Distribution of \ 2 f° r improved sta- 
bility assumption (equation 9.6). Solid 
line is the best fit of \ 2 distribution func- 
tion with 337.6 degrees of freedom. The 
goodness of this fit is characterized by 
X 2 /ndf = 92.33/92. 



Figure 9.7: Plots showing distributions of \ 2 of the difference of simulated and real 
azimuth angle distributions for every half hour with and without correction. The 
expected number of degrees of freedom is 337.5. Both histograms have the same 
number of entries (2073). 
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90 bins, and azimuth one into 360. A word of precaution is appropriate: the real 
and generated distributions are not independent, one is produced from the other 
using swapping technique. Due to implementation of the time swapping method (8 
hour integration time window) and because distributions are accumulated every 30 
minutes, the number of degrees of freedom in the y 2 test should be lower than one 
would expect from the usual histogram bin considerations by 1 for every 16 bins or 
84.25 and 337.5 for 90 and 360 bin histograms respectively. This is supported by 
Monte Carlo simulations presented on figure (9.5) where the values of x 2 resulted 
from the test are histogramed. 

Figure (9.6) demonstrates zenith angle distribution comparisons with and with- 
out the correction. From these plots it is clearly seen how a bias of the standard 
method is removed. A similar study was carried out for azimuth angle distribu- 
tion. Figure (9.7) demonstrates the results of the study and shows that azimuth 
variation, if any, is on a much smaller scale than the zenith one. Note, that no 
azimuthal correction was applied, therefore, a slight change in the azimuth dis- 
tribution test is due to possible correlation between zenith and azimuth angles of 
the detected events. Thus, the improved stability assumption (equation 9.6) is 
considered to be valid for the 8 hour duration of timing integration. 

9.7 Known Anisotropics. 

It was, so far, conceded that no anisotropy on the sky is present. This, together 
with the stability assumption had lead to the equation (9.1). In fact, if there are 
known sources on the sky, then the number of registered events is given by: 

dN(x, t) = S(x, t) ■ G(x) ■ R(t) dxdt 

where S(x, t) describes the strength of the sources as function of local coordi- 
nates and time. (This is how anisotropy can mimic detector instability.) 

There are two known sources, actually sinks, on the sky: the Sun and the Moon 
[57]. Because the anisotropy function S(x,t) produced by them is not known, the 
anisotropy is handled by vetoing the ±5° degree regions around the objects. This 
entails treating them as part of the source bin during Monte Carlo integration 
and disregarding generated events in the background estimation 7V& if they fall 
within the veto region. In general, known small scale anisotropies must be vetoed 
as described, known large scale ones have to be incorporated into the stability 
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assumption. These become part of the null hypothesis. 
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Chapter 10 

Alternative Hypothesis. 



The construction of the alternative hypothesis, or more precisely, of critical regions 
is based on the knowledge about the Galactic gamma ray emission available to date 
[9, 7]. Among all experiments operated at different energies the EGRET telescope 
is the one of closest energy range to Milagro and made a convincing detection 
of emission. The data made available by the EGRET collaboration [59, 7] is the 
foundation of the alternative hypothesis. Because of the change of the primary 
emission mechanism between MeV and GeV energy regions from bremsstrahlung to 
7T° decay [14], only highest EGRET energy data points are used in the formulation. 
Figure 10.1(a) plots the longitude dependence of the EGRET emission flux in the 
±2° latitude band along the Galactic equator at energies 10-30 GeV. Figure 10.1(b) 
shows the Milagro Galactic disk exposure as a function of Galactic longitude. The 
expected significance which is roughly proportional (equation (8.1)) to the product 
of the signal flux and the square root of the exposure is presented on figure 10.1(c). 

Thus, the region of Galactic disk with I G (20°, 100°) is the region where highest 
significance is expected provided the emission model of EGRET is valid at TeV 
energies. Therefore, the region of I G (20°, 100°) is selected as Milagro inner 
Galaxy. Due to the fact that the inverse Compton mechanisms may result in a 
wider latitude distribution than that of the gas column density, the region with 
\b\ < 5° is also considered. The two regions with I G (140°, 220°), \b\ < 2° and 
|6| < 5° are also selected as Milagro outer Galaxy. These parts of the disk are 
selected due to their high exposure. These regions are analyzed separately, however 
|6| < 2° and \b\ < 5° ones are not independent for the same longitude range. The 
regions were selected before the analysis was performed. 
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(a) Total flux at 10-30 GeV in |6| < 2° disk. 
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(b) Milagro exposure of \b\ < 2° disk. 
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(c) The expected relative significance. 



Figure 10.1: Plots showing the construction of the critical region for the alternative 
hypothesis based on the EGRET data. 
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Chapter 11 
Results. 



11.1 Results of the Measurements. 

The results presented below are obtained after the analysis of data collected be- 
tween July 19, 2000 and September 10, 2001 by the Milagro detector. The starting 
date corresponds to the commissioning of the hadron rejection algorithm, the end- 
ing date is when the experiment was turned off for a major scheduled maintenance. 

To ensure high quality of the data, the periods with abnormal event rate with 
and without the hadron rejection cut were eliminated from the analysis. These 
may indicate possible detector problems. Parts of the data were also discarded 
if the Milagro log-book showed entries of some hardware/software problems or 
indications of unstable operation. Periods of abnormal zenith angle distribution 
of events were also discarded. Other cuts applied include Nf it > 20 to ensure 
good angular resolution and the zenith angle cut 0° < z < 50° due to stability 
considerations (appendix C). 

The null hypothesis tested consists of: 

• There is no emission from the Galactic plane except for background. 

• All of the detected events are due to background and it is isotropic except 
for ±5° region around the Sun and the Moon. (Events from these regions 
are vetoed.) 

• Detected events satisfy the stability equation (9.6). 
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• The null hypothesis also includes certain aspects of shower development 
and detector operation important for the detection and reconstruction of 
events. These include reconstructions algorithms and detector calibration 
procedures. 

Because the emission from the Galactic plane is under investigation, events 
arriving from the ±7° band along the Galactic equator are excluded from the 
background generation (see section 9.2). Four regions of the disk (selected a priory) 
were studied. The Milagro inner Galaxy with Galactic longitude I e (20°, 100°), 
\b\ < 2° and \b\ < 5°; and Milagro outer Galaxy with I E (140°, 220°), \b\ < 2° and 
|6| < 5°. The results are given in the table 11.1. 





Inner Galaxy 


b < 2° 


b <5° 


N s 


(43446.6 ± 6.6) • 10 3 


(10823.8 ± 1.0) • 10 4 


N b 


(43426.7 ±2.7) • 10 3 


(10819.7 ±0.5) • 10 4 




(±4.57 ± 1.64) • 10" 4 


(±3.77 ± 1.08) • 10" 4 


significance 


±2.8 


±3.5 




Outer Galaxy 


b < 2° 


b <5° 


N s 


(46409.1 ±6.8) • 10 3 


(11588.6 ± 1.1) • 10 4 


N b 


(46413.0 ± 3.0) • 10 3 


(11590.1 ±0.6) • 10 4 


NjN b 


(-0.82 ± 1.6) • 10~ 4 


(-1.28 ± 1.06) • 10" 4 


significance 


-0.5 


-1.2 



Table 11.1: Summary of the results for the analysis of the Milagro inner and outer 
Galaxy (JV 7 = N s - N b ). 

As seen from the table 11.1 the results of the analysis for the inner Galaxy 
indicate violation of the null hypothesis. The null hypothesis is rejected with the 
significance of 3.5 or the probability that it is rejected by chance fluctuations is 
only 2.3 • 10~ 4 (see table 8.1). A detection of a source is usually claimed if the 
chance probability is less than 0.1%. 

In order to be able to interpret the result as a presence of excess of events in 
the direction of the Galactic plane, other constituents of the null hypothesis have 
to be examined. The aspects of detector operation and event reconstruction are 
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unlikely to produce localized structure on the sky, they should rather show up as 
global changes such as in total event rate, in zenith angle distribution of detected 
events. These were carefully monitored. A violation of the stability assumption 
may, however, "produce" large scale structures on the sky (and vise versa). 

Isotropy of the background is another ingredient. Because of the random bend- 
ing of cosmic rays in the galactic magnetic fields which are responsible for the 
isotropy of the background, small deviations from isotropy are expected to be 
smooth and therefore reveal themselves as large scale structures too. According 
to indications in the literature from underground muon experiments [60] a slight 
anisotropy may be present which is parametrized by the authors as: 

R(a) = 1 + r cos(« — a ) 

where r = (5.6 ± 1.9) ■ 10~ 4 , = 8.0° ± 19.1°. To study the systematic error 
caused by a possible neglect of the effect, a computer simulation was employed. In 
this simulation a data set equivalent to the 18 times the real one was generated. 
The anisotropy was simulated at the specified level but standard analysis was 
applied. The results of this simulation are presented in the table 11.2. 





\b\ < 2° 


b <5° 


Inner Galaxy 


(-3.7 ±3.8) • 10- 5 


(-4.0 ±2.5) • 10~ 5 


Outer Galaxy 


(-5.4 ±3.7) • 10- 5 


(-1.7 ±2.4) • 10" 5 



Table 11.2: Results of the systematic error study for N 1 /N h . 

Despite the computational error bars that are only about 4.3 times smaller than 
the ones in the real data (limited by the computer time), this is enough to show 
that this anisotropy is also unlikely to be responsible for the observed result. It 
should be noted that the discussed muon anisotropy need not be identical to that 
of the cosmic ray air showers in the Milagro energy range. However, the result 
may be interpreted as a systematic correction if the anisotropy is established to 
exist [61]. Within limits of the calculation, the correction is not inconsistent with 
being zero. 

In the process of background estimation using the time swapping method, the 
background was generated and saved for the entire sky, not only for the regions of 
interest. This may serve to provide additional information. 
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Figure 11.1: Significance map of the sky in Galactic coordinates. Emission from 
the inner Galactic region is seen. Grid lines are plotted every 30° in longitude and 
latitude. The Galactic center, not visible by Milagro, is in the middle of the map. 
Galactic longitude increases to the left. 
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Figure 11.2: The distribution of significances of individual cells of figure 11.1. The 
best fit of a Gaussian distribution is also shown. It agrees well with the standard 
normal distribution (zero mean and unit variance). 
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(a) Latitude profile of N 7 /Nb for inner Galaxy. The value of x 2 with respect to 
is 72.10 with 83 bins. 
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(b) Latitude profile of Nj/Nb for outer Galaxy. The value of x 2 with respect to 
is 116.9 with 84 bins. 



Figure 11.3: Latitude profiles of Nj/N b for the inner and outer Galaxy. 



93 



Figure 11.1 shows the 2 dimensional map of the significance in the Galactic 
coordinates with coarse bins appropriate for a large scale structure. An enhance- 
ment in the region of inner Galaxy is seen. Another feature visible on the plot is a 
valley at I = —150°, b = +60°. It may be interpreted as an indication of a possible 
small anisotropy of the cosmic rays. Haverah Park [62] and Yakutsk [63] experi- 
ments have found a deficit of flux from the Galactic North in the energy range of 
10 15 — 10 19 eV, however, the existence of such an anisotropy in the Milagro energy 
range is unknown. Figure 11.2 presents the distribution of significances of individ- 
ual cells of figure 11.1. The distribution is expected to be the standard normal one 
(see discussion following equation (8.1)) with which it agrees well. Therefore, the 
observed excesses and deficits are statistically allowed. 

A signature of the galactic nature of the emission is presence of a ridge along 
the galactic equator. Therefore, another check is provided by a lateral distribu- 
tion plot of N y /N b (figure 11.3). The plot for the inner Galaxy (figure 11.3(a)) 
shows the excess in the region around the equator (b = 0) as expected and is flat 
outside suggesting that it is the disk which produces the excess. The plot for the 
outer Galaxy (figure 11.3(b)) does not show presence of excess in the plane region, 
however, indicates a possible existence of a large scale anisotropy. The exact in- 
terpretation of the values of x 2 given in the captions to the figures is, however, 
difficult. Calculation of values and errors for the points plotted outside of the 
excluded region suffers from the dependence problem such as illustrated on figure 
9.1. Also, the points in adjacent bins are positively correlated, therefore, number 
of degrees of freedom is effectively lower than what is expected from the number 
of bins considerations. However, if it is assumed that anisotropy is present and 
is modeled by a straight line in the inner Galactic region and by a parabola in 
the outer one, then it is possible to estimate its effect by fitting a straight line 
in the region of 5° < |6| < 60° for the case of the inner Galaxy and quadratic 
polynomial in the region of 5° < |6| < 40° for the case of the outer Galaxy and 
interpolating it into the plane \b\ < 5°. This results in the estimated contribution 
of (+0.54 ±0.34) • 10~ 4 (inner Galaxy) and (-0.48 ±0.65) • 10~ 4 (outer Galaxy) to 
the measured ratios N 1 /N cr . These values have to be subtracted from the results 
for the ratios of table 11.1 as systematic corrections if the anisotropy is known to 
be of the form discussed [61]. They also are not inconsistent with being zero. 

It is therefore concluded that the outer Galaxy does not show excess and the 
measurements may be used to set upper limits for the emission. The results for 
the inner Galaxy, if interpreted as gamma ray emission, may be used to calculate 
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the gamma ray flux. 



11.2 Interpretation of the Results. 

In what follows it is assumed that the emission from the Galactic plane is known 
to exist and that the observed excess in the galactic plane is due to gamma rays. 
With these assumptions in mind, the results of the measurements can be used to 
establish some properties of the emission. 

Assuming simple power law emission spectrum independent of Galactic coor- 
dinates for Galactic gamma rays in the Milagro energy range or, equivalently, 
assuming gamma ray signal function to be in the form of equation (8.4) and 
given power law spectrum of cosmic-ray flux at the Earth (F cr (> ITeV) — 1.2 • 
1CT 5 (cm~ 2 s~ 1 str~ 1 ) [64], a cr = 2.7) it is possible to use results of the table 11.1 to 
relate integral flux of gamma rays above 1 TeV (F 7 (> ITeV)) with their assumed 
spectral index a 7 (equation 8.5). The dimensionless coefficient r](a y ,a cr ) needed 
for this conversion is plotted on figure 11.4 for a cr = 2.7. The value of this coeffi- 
cient reflects efficiency of the cosmic ray rejection algorithm and for a 1 = a cr has 
a direct interpretation of a ratio of effective areas after the cut. The integral flux 
of gamma rays above 1 TeV obtained is plotted on figure 11.5 as solid lines. 





\b\ < 2° 


\b\ < 5° 


Inner Galaxy 


(1.82 ±0.25) • 10-* 


(1.18 ±0.14) • 10~ s 


Outer Galaxy 


(1.90 ±0.81) • 10~ 7 


(3.45 ±0.60) • 10~ 7 



Table 11.3: Average integral flux of photons (cm~ 2 sr^ 1 s^ 1 ) in the 10-30 GeV 
energy range measured by EGRET in the corresponding regions of the Galactic 
plane. Isotropic diffuse component has been subtracted [7, 59]. 

On the other hand, given the measurements of EGRET (table 11.3) and assum- 
ing that power law index a 7 is constant above 10 GeV it is possible to extrapolate 
the EGRET measurement in the highest energy bin 10-30 GeV to the 1 TeV re- 
gion. This is shown on figure 11.5 as dashed lines. Both EGRET and Milagro 
measurements are presented together with the corresponding error corridors. The 
intersection of the corridors provides the spectral index and integral flux that is 
consistent with the EGRET and Milagro data. In the case of Milagro upper limit, 
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Figure 11.4: Dimensionless coefficient ?7(a: 7 ,a: cr ) plotted as function of a 7 (a cr = 
2.7) for the inner Galaxy with |6| < 5°. Other regions of the Galaxy considered have 
very similar exposure and therefore the presented coefficient is almost identical. 
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(a) Inner Galaxy, |6| < 2°. 



(b) Inner Galaxy, \b\ < 5°. 
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(c) Outer Galaxy, |6| < 2°. 



(d) Outer Galaxy, |6| < 5°. 



Figure 11.5: Integral flux of photons in different parts of Galactic disk with energies 
above 1 TeV or corresponding 99.9% upper limit for assumed differential spectral 
index a 7 are shown as solid lines. The dashed lines represent the extrapolated 
EGRET measurements in the 10-30 GeV energy range and the same filed of Galaxy. 
The error corridors are shown as lines with dots. 
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the limit itself is presented. This comparison leads to estimation of the spectral 
index a 7 = 2.59 ± 0.07 and flux F(> ITeV) = (9.5 ± 2.0) • 10" 10 cm' 1 sr' 1 s' 1 
in the region of inner Galaxy and to a lower limit of a 7 > 2.49 and flux F(> 
ITeV) < 4.5 • 10~ 10 cm' 1 sr' 1 s" 1 in the outer Galaxy. (Regions with |6| < 5° are 
used for final calculations because of greater sensitivity of Milagro.) These find- 
ings may be used to exclude models which predict a strong enhancement of the 
diffuse flux compared to conventional mechanisms. Model predictions attempting 
to explain the excess flux in the GeV region measured by EGRET by assuming 
an increased inverse Compton component [13] or by harder proton spectrum [1] 
or by considering contribution from unresolved supernova remnants [2] are given 
in literature for different ranges in Galactic latitude and longitude and can not be 
directly compared with the presented results. However, it is possible to find spec- 
tral indexes in the energy range above 1 GeV as measured by EGRET averaged 
over the Milagro inner and outer Galaxies (table 11.4) and compare them with the 
above results. While the EGRET spectral index does not contradict to the Milagro 
upper limit for the region of outer Galaxy, measurements of the index in the inner 
Galaxy exclude the possibility of constant spectral index and require softening of 
the spectrum between EGRET and Milagro energy regions. This is also illustrated 
on figures 11.6 for the regions of inner and outer Galaxy with \b\ < 5°. 





\b\ < 2° 


\b\ < 5° 


Inner Galaxy 


2.346 ± 0.028 


2.365 ±0.016 


Outer Galaxy 


2.560 ±0.054 


2.487 ±0.059 



Table 11.4: Spectral index a 7 derived from the EGRET data in the 1-30 GeV 
energy range in the corresponding regions of the Galactic plane. (The data was 
made available by [59].) 
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(a) Inner Galaxy. 



(b) Outer Galaxy. 



Figure 11.6: Integral flux of photons from \b\ < 5° region multiplied by photon 
energy E is plotted. EGRET data points in the energy range of 1-30 GeV are also 
shown. The dashed line is the spectrum derived from the EGRET data above 1 
GeV. 
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Chapter 12 
Conclusions. 



A search for diffuse gamma ray emission from the Galactic plane has been per- 
formed using data collected by the Milagro detector between July 19, 2000 and 
September 10, 2001. Four regions of the disk were studied. An excess has been 
observed from the regions of the Milagro inner Galaxy defined by I G (20°, 100°), 
\b\ <2°and|6| < 5° with the significance of 2.6- 10~ 3 and 2.3- 10~ 4 respectively. The 
emission from the region of the Milagro outer Galaxy defined by I G (140°, 220°), 
|6| < 2° and |6| < 5° is not inconsistent with being that of background only. 

The results are interpreted assuming a simple power law energy spectrum. 
Given the EGRET measurements, the integral gamma ray flux with energies above 
1 TeV has been obtained, yielding F(> ITeV) = (9.5 ± 2.0) • 10~ 10 cmr 1 sr^ 1 s^ 1 
with spectral index a 7 = 2.59 ± 0.07 for the region of inner Galaxy. The 99.9% 
upper limit for the diffuse emission in the region of outer Galaxy is set at F(> 
ITeV) < 4.5 • 10~ 10 cvrT x sr^ 1 s^ 1 assuming a differential spectral index of a 1 = 
2.49. These observations are compatible with earlier upper limits obtained by the 
Whipple observatory at 500 GeV [8]. The upper limit for the outer Galaxy is 
also consistent with the extrapolation of EGRET measurements between 1 and 30 
GeV. Extrapolation of the EGRET measurements between 1 and 30 GeV for the 
region of inner Galaxy using constant power law spectral index is incompatible 
with the Milagro data. This indicates softening of the spectrum at energy between 
10 GeV and 1 TeV. These observations may be used to constrain models of Galactic 
emission [1, 2, 13]. 

During preparation of this manuscript, the Milagro detector continued to oper- 
ate and collected more data. The analysis of these data will confirm or refute the 



100 



observations. Regardless of the outcome, the combined data set will provide for 
more sensitive measurement of the diffuse gamma ray emission at energies near 1 
Tev. The Milagro detector is being continuously upgraded both on software and 
hardware level. New event reconstruction algorithms may become available in the 
near future. These, however, will have no affect on the data collected so far. In 
terms of hardware, Milagro is being supplimented by an outrigger array, which 
will enable better hadron rejection as well as energy determination and angular 
resolution. As more and better quality data become available, it will be possible to 
study the question of cosmic ray anisotropy extensively using the Milagro data it- 
self. This presents an interest not only in the context of the search for the Galactic 
gamma ray emission, but is valuable by itself as the detection of Compton-Getting 
effect may become possible [65]. 
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Appendix A 



Arguments pertaining to equation 
8.1. 



The Poisson distribution of an integer random variable n is given by: 

P» = £e-" 

where fj, is a parameter of distribution. One can show that both average and 
dispersion of n are equal to /i. It is intended to show that for large (j,, the Poisson 
distribution approaches Gaussian: 

I _(fzfo)! 



27TCT 

with average xq = fj, and dispersion cr 2 = //. 

For the case of large \x and n being not too different from it, the Stirling formula 
for the factorial can be applied: 



n\ pa v / 27m n n e n 



with relative error of order of (eian — 1) ; leading to 



\j2ixn \n 
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The exponent y(n) can be expanded in Taylor series around \i: 



v(p) = 


— In y/27V fJ, 


y'(») = 


1 


y"{v) = 




y"'{n) = 









In this Taylor series (recalling that \x is large) the second term can be neglected 




y(n) = - In y^ - -L(„ - „) - -L h + i-J (n - /i) 2 + 



compared to the leading if |n — //| <C 2/x In y/2ir[/,, third if |n — /j| <C y 2/x In y/2ir/i, 
fourth if \n— fi\ <C (6/i 2 In ^/2ixji) 1 l z , all higher terms if \n— fx\ <C In y/2nfj,y/ k . 
It is seen that the first condition is least restrictive one, second, however, is the 
most limiting. Thus, keeping only two leading terms in \n — fx\ will result: 



y(n) « - In y^ - J- ( l + J- ) ( n - ^) 



with restriction \n — fj,\ <C (6/i 2 In v / 27r / u) 1//3 . This means that ra has to be not 
too far off \x and when /x is large the use of the Sterling formula is justified. The 
Poisson distribution is approximated by Gaussian in the central region: 

P^(n) = ^-e'^ -^=e~^^, \n - n\ < (6fi 2 In y^) 173 , n > 3/2 
n! yz7r/x v 

Examining the equation 8.1 if null hypothesis is true, it follows from the above, 
that the numerator is distributed under the Gauss law with zero mean if both Ni 
and N 2 are large and are not too different from = Xt i;2 . By construction, if 
null hypothesis is true and N^ 2 ~ /ii,2, the dispersion of S is approximately unity, 
therefore, the statistic S has normal distribution in — So < S < So region around 
zero. The bound So is defined by the smallest of the iV's and is on the order of 
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(a) Null hypothesis with parameters \xi = 
5000, fii = a^2, a = 0.1 was used. Ac- 
cording to equation (A.l), So <C 8. 



(b) Null hypothesis with parameters jj.2 = 
15 • 10 6 , jUi = a/U2, a = 0.1 was used. 
According to equation (A.l), So <C 39. 



Figure A.l: Distributions of statistic S (equation 8.1) obtained in two runs of 
Monte Carlo simulations. Best fit to Gauss distribution is made, parameters of 
the fit are presented. The two curves should agree in the So-neighborhood around 
zero. Number of entries in each run (about 10 9 ) was chosen to provide reasonable 
accuracy in the region plotted. 



A Monte Carlo simulation helps to verify the statement. Figure A.l shows the 
results of such a simulation. 




(A.l) 
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Appendix B 



Solution of Background 
Equations. 



The problem of exclusion of the source region from the background estimation had 
led to the problem of solution of the integral equations: 



with respect to G(x) and R(t). N out (x) and R ou t{t) represent the input data 
collected from the off-source region, <f>(x, t) is the function which defines its bounds: 
it is one if (x, t) points into the off-source region and zero otherwise. It is first 
noted that both G(x) and R(t) enter into equations only as a product G(x) ■ R(t), 
therefore, normalization of either of them does not make any difference as long 
as the product is preserved. Also, if there is a point xq in the local coordinates 
which is never exposed into the off-source region, that is (f>(xo,t) = 0, Vt, then 
N out (x ) = and the first equation becomes: 



leading to G(xq) being undefined. Such regions, if any, are marked with a nega- 
tive value, on-source events with local coordinates from these regions are discarded 
as having no corresponding background estimate. In such cases, measurement is 
impossible within the environs of assumptions. 

The system of equations is solved numerically. The equations are discretized 




G(x) f <p(x,t') -R{t') dt' 
R(t) J 4>{x' ,t) ■ G(x ) dx' 



= G(x ) ■ 
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® Declination 



Hour Angle 



Figure B.l: To the formulation of backgroud equations. 
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Figure B.2: To the derivation of recursion relations. 
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on a fine grid corresponding to 0.1° in x and t, the integrals become sums, x 
and t become indexes. Since most often one is interested in celestial objects, the 
off-source region is defined in declination and right ascension (5, a) and due to 
discretization presents the bound indexes a± and 0,2-, ol\ < ol 2 for every declination 
index 5. Therefore, for this problem, it is natural to consider local coordinates x 
as declination and hour angle (S, h) rather than zenith and azimuth. The situation 
is illustrated on figure B.l. The shaded area is the off-source region bounded by 
(f>(x,t) in its discrete form. The region of interest, the on-source region, is defined 
by some other conditions which are irrelevant for the background equations as long 
as both regions do not overlap. 

Thus, the excluded region, defined by interval (ai,a 2 ) is translated to local 
coordinates h and t as region inside intervals (hi, h 2 ) and (ti,t 2 ), where 

hi = t — a 2 mod 360° 
h 2 = hi + a 2 — aci 

ti = h + a 2 mod 24 hours 
h = h + a 2 — a>i 

Because of these relations it is convenient to quantize all quantities in the same 
fashion. This, however, is not required. Here, it is adopted that a\ is the index of 
the first excluded bin, a 2 is the index of the last excluded bin. So will be hi,ti and 
h 2 ,t 2 . The integration is performed using recursion formulas which are illustrated 
on the p(t) = J <f)(h',t)G(h')dh' example. After integral is replaced by sum, it 
becomes: 

hi-l N-l 

p{t) = £ G(h) + £ G(h) 

h=0 h=h 2 +l 

where N is the number of cells from to iV — 1 created by the discretization 
of the hour angle between 0° and 360°. The limits of the summation are included, 
that is 

b 

£ G(h) = G(a) +G(a+1) + --- + G(b-1) + G(b) 

h=a 

Then, it is seen (figure B.2) that the integral at some other time t' is given by 



107 



h[-l N-l 

p(t') = ^ G(h) + £ = p(t) + Ax - A 2 

/i=0 h=h' 2 +l 



where 




Ai=+e!£ ^(ft), K>hi 
A 1 = -T,lr 1 G(h), h[<hi 



Ai = 0, /ii = /ii 



A 2 = +E^ +1 G(/ i ), h 2 >/i 2 
A2 = -E^ +1 G(/i), ^</i 2 
A 2 = 0, h' 2 = h 2 



This scheme is applied for every declination. Integration over time is performed 
in the similar fashion. It is important to note, that case of a± = a 2 corresponds 
to exclusion of one bin. Absence of excluded region for a declination is encoded as 
0(2 — Oil — 1) which is the only case with a 2 < a\ allowed. The system is solved 
by iterations where R out {t) is used as first approximation to R(t). Then, the first 
equation is solved with respect to G(x) which is used in the second equation to 
update R(t). The process continues until desired accuracy is reached. In practice, 
13-digit precision is achieved after about 15 iterations. 
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Appendix C 
Zenith Correction. 



The zenith correction function A(z) describes deviation of the current shape of the 
zenith distribution from the average one over the same UT day. It is obtained in 
the iteration process in the form of polynomial of 5-th degree starting with some 
zero-th approximation. The differences between the two normalized histograms 
are accumulated according to the sign of amplitude 9 of the fit of a difference into 
9A(z). When all differences according to their corresponding sing of amplitudes 
are accrued, coefficients of the A(z) are adjusted to provide the best fit. This 
is the updated function A(z). The procedure is repeated until coefficients of the 
polynomial stop changing. Then, A(z) is rescaled so that its maximum is equal to 
one. This is the final form of A(z). The K(z) function is obtained from the A(z) 
by 

K(z) - A ^ 



G(z) 



where G(z) is the average zenith distribution. Coefficients of the correction 
function A(z) are given in the table C.l with its shape defined by: 

A(z) = A(y 5 + ey A + dy 3 + cy 2 + by + a) 
where y = z/90, 0° < z < 50°. 

The example of the A(z) curve is presented on figure C.l The field of view of 
the detector was artificially made limited to 0° < z < 50° because zenith angle 
distribution variations outside of this range could not be characterized by such a 
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b 




A 






(MJD) 


(io 3 - 6 ) 


(10- 3 ) 


(IO- 4 ) 


(io- 1 ) 


(io"- 2 ) 


(10°) 


1745 














I 


+203.0670 


-174.6125 


+10721.45 


-16.27103 


-5.631589 


-122.02 


1813 














1 


-27.00866 


-19.20310 


+291.3698 


+4.232470 


-129.1226 


-698.0 


1876 














1 


+56.50721 


-32.53573 


+1355.857 


+ 1.621818 


-108.8645 


-551.0 


1937 














I 


+124.5009 


-12.71605 


+ 171.8406 


+3.488231 


-116.1735 


-1300.0 


2000 














I 


-4.604483 


-10.42552 


+4.944296 


+4.006497 


-121.4843 


-1260.0 


2069 














I 


-32.60869 


-8.961481 


-297.4525 


+5.175355 


-133.4420 


-1100.0 


2117 














I 


-27.23581 


-11.71815 


+27.33091 


+4.266234 


-126.5341 


-1110.0 


2163 















Table C.l: Coefficients of the zenith correction function A(z) derived for data 
which passes the hadron rejection cut X2 > 2.5. MJD, the Modified Julian Day, 
identifies the validity time period. For instance, MJD 1745 corresponds to July 
19, 2000; MJD 2163 corresponds to September 10, 2001 [58]. 
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Figure C.l: Correction function A(z). Coefficients corresponding to period MJD 
1745 - MJD 1813 were used. 
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